Changeset 15134


Ignore:
Timestamp:
2021-07-22T18:26:16+02:00 (3 months ago)
Author:
dbruciaferri
Message:

Cleaning and reorganising the code. It replicates JMMP-AMM7 domain_cfg_MEs_L51_r24-07_opt_v2.nc.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/mes.F90

    r15131 r15134  
    44   !! Ocean initialization : Multiple Enveloped s coordinate (MES) 
    55   !!============================================================================== 
    6    !!  NEMO      3.6  ! 2017-05  (D. Bruciaferri)    
     6   !!  NEMO      3.6   ! 2017-05  D. Bruciaferri 
     7   !!  NEMO      4.0.4 ! 2021-07  D. Bruciaferri   
    78   !!---------------------------------------------------------------------- 
    8  
    9    !!---------------------------------------------------------------------- 
    10    !!   dom_mes          : define the MES vertical coordinate system 
    11    !!       fs_ma96      : tanh stretch function Madec 1996 
    12    !!       fs_sh94      : Song and Haidvogel 1994 stretch function 
    13    !!--------------------------------------------------------------------- 
     9   ! 
    1410   USE oce               ! ocean variables 
    1511   USE dom_oce           ! ocean domain 
     
    3531   INTEGER            :: nn_strt(max_nn_env)     ! Array specifing the stretching function for each envelope: 
    3632                                                 ! Madec 1996 (0), Song and Haidvogel 1994 (1)  
    37    REAL(wp)           :: rn_ebot_max(max_nn_env) ! depth of envelopes in one point (for checking 1D transform.) 
    3833   REAL(wp)           :: max_dep_env(max_nn_env) ! global maximum depth of envelopes 
    3934   REAL(wp)           :: min_dep_env(max_nn_env) ! global minimum depth of envelopes 
     
    6762      !! 
    6863      !! ** Method  :   s-coordinate defined respect to multiple 
    69       !!                externally defined envelope 
    70       !! 
    71       !!      The depth of model levels is defined as the product of an 
    72       !!      analytical function by the externally defined envelopes. 
     64      !!                externally defined envelope as detailed in 
     65      !! 
     66      !!                Bruciaferri, Shapiro, Wobus, 2018. Oce. Dyn.  
     67      !!                https://doi.org/10.1007/s10236-018-1189-x 
     68      !! 
     69      !!                Three options for stretching are given: 
    7370      !!  
    74       !!      - Read envelopes (in meters) at t-point and compute the 
    75       !!        envelopes at u-, v-, and f-points. 
    76       !!            hbatu = mi( hbatt ) 
    77       !!            hbatv = mj( hbatt ) 
    78       !!            hbatf = mi( mj( hbatt ) ) 
    79       !! 
    80       !!      - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical 
    81       !!        function and its finite difference derivative. 
    82       !!            z_gsigt(k) = fssig (k    ) 
    83       !!            z_gsigw(k) = fssig (k-0.5) 
    84       !!            z_esigt(k) = z_gsigw(k+1) - z_gsigw(k) 
    85       !!            z_esigw(k) = z_gsigt(k+1) - z_gsigt(k) 
    86       !! 
    87       !!      Three options for stretching are given: 
    88       !!  
    89       !!          *) nn_strt = 0 : Madec et al 1996 cosh/tanh function 
    90       !!          *) nn_strt = 1 : Song and Haidvogel 1994 sinh/tanh function   
    91       !!          *) nn_strt = 2 : Siddorn and Furner gamma function 
    92       !! 
    93       !!---------------------------------------------------------------------- 
    94       !!        SKETCH of the GEOMETRY OF THE S-S-Z COORDINATE SYSTEM 
     71      !!                *) nn_strt = 0 : Madec et al 1996 cosh/tanh function 
     72      !!                *) nn_strt = 1 : Song and Haidvogel 1994 sinh/tanh function   
     73      !!                *) nn_strt = 2 : Siddorn and Furner gamma function 
     74      !! 
     75      !!---------------------------------------------------------------------- 
     76      !!        SKETCH of the GEOMETRY OF A MEs-COORDINATE SYSTEM 
    9577      !! 
    9678      !!        Lines represent W-levels: 
     
    138120 
    139121      ! 
    140       INTEGER                             ::   ji, jj, jk, jl, je       ! dummy loop argument 
    141       INTEGER                             ::   num_s, s_1st, ind        ! for loops over envelopes 
    142       INTEGER                             ::   num_s_up, num_s_dw       ! for loops over envelopes 
    143       REAL(wp)                            ::   D1s_up, D1s_dw           ! for loops over envelopes 
    144       INTEGER                             ::   ibcbeg, ibcend           ! boundary condition for cubic splines 
    145       INTEGER, PARAMETER                  ::   n = 2                    ! number of considered points for each interval 
    146       REAL(wp)                            ::   c(4,n)                   ! matrix for cubic splines coefficents 
    147       REAL(wp)                            ::   d(2,2)                   ! matrix for depth values and depth derivative  
    148                                                                         ! at intervals' boundaries 
    149       REAL(wp)                            ::   tau(n)                   ! abscissas of intervals' boundaries 
    150       REAL(wp)                            ::   coefa, coefb             ! for loops over envelopes 
    151       REAL(wp)                            ::   alpha, env_hc            ! for loops over envelopes 
    152       REAL(wp)                            ::   dep_up, dep_dw           ! for loops over envelopes 
    153       REAL(wp)                            ::   dep0, dep1, dep2, dep3   ! for loops over envelopes  
    154                                                                         ! for cubic splines 
    155       REAL(wp)                            ::   zcoeft, zcoefw, sw, st   ! temporary scalars 
    156       REAL(wp)                            ::   max_env_up, min_env_up   ! to identify geopotential levels 
    157       REAL(wp)                            ::   max_env_dw, min_env_dw   ! to identify geopotential levels 
    158       REAL(wp), POINTER, DIMENSION(:)     ::   z_gsigw1, z_gsigt1       ! MES-coordinates analytical trans.  
    159                                                                         ! (0.<= s <= 1.) for T- and W-points 
    160       REAL(wp), POINTER, DIMENSION(:)     ::   z_esigw1, z_esigt1       ! MES analytical scale factors   
    161                                                                         ! of T- and W-points 
    162       REAL(wp), POINTER, DIMENSION(:,:  ) :: env_up, env_dw             ! for loops over envelopes 
    163       REAL(wp), POINTER, DIMENSION(:,:  ) :: env0, env1, env2, env3     ! for loops over envelopes 
    164                                                                         ! for cubic splines 
     122      INTEGER                               :: ji, jj, jk, jl, je       ! dummy loop argument 
     123      INTEGER                               :: eii, eij, irnk           ! dummy loop argument 
     124      INTEGER                               :: num_s, s_1st, ind        ! for loops over envelopes 
     125      INTEGER                               :: num_s_up, num_s_dw       ! for loops over envelopes 
     126      REAL(wp)                              :: D1s_up, D1s_dw           ! for loops over envelopes 
     127      INTEGER                               :: ibcbeg, ibcend           ! boundary condition for cubic splines 
     128      INTEGER, PARAMETER                    :: n = 2                    ! number of considered points for each  
     129      !                                                                 ! interval 
     130      REAL(wp)                              :: c(4,n)                   ! matrix for cubic splines coefficents 
     131      REAL(wp)                              :: d(2,2)                   ! matrix for depth values and depth  
     132      !                                                                 ! derivative at intervals' boundaries 
     133      REAL(wp)                              :: tau(n)                   ! abscissas of intervals' boundaries 
     134      REAL(wp)                              :: coefa, coefb             ! for loops over envelopes 
     135      REAL(wp)                              :: alpha, env_hc            ! for loops over envelopes 
     136      REAL(wp)                              :: zcoeft, zcoefw           ! temporary scalars 
     137      REAL(wp)                              :: max_env_up, min_env_up   ! to identify geopotential levels 
     138      REAL(wp)                              :: max_env_dw, min_env_dw   ! to identify geopotential levels 
     139      REAL(wp)                              :: mindep 
     140      ! 
     141      REAL(wp), DIMENSION(max_nn_env)       :: rn_ebot_max 
     142      REAL(wp), DIMENSION(jpk)              :: z_gsigw1, z_gsigt1       ! 1D arrays for debugging 
     143      REAL(wp), DIMENSION(jpk)              :: gdepw1, gdept1           ! 1D arrays for debugging 
     144      REAL(wp), DIMENSION(jpk)              :: e3t1, e3w1               ! 1D arrays for debugging 
     145      ! 
     146      REAL(wp), DIMENSION(jpi,jpj)          :: env_up, env_dw           ! for loops over envelopes 
     147      REAL(wp), DIMENSION(jpi,jpj)          :: env0, env1, env2, env3   ! for loops over envelopes 
     148      !                                                                   ! for cubic splines 
     149      REAL(wp), DIMENSION(jpi,jpj)          :: pmsk                     ! working array 
     150      ! 
     151      REAL(wp), DIMENSION(jpi,jpj,jpk)      :: z_gsigw3, z_gsigt3       ! 3D arrays for debugging 
    165152      !!---------------------------------------------------------------------- 
    166153      ! 
     
    169156      CALL mes_init 
    170157      ! 
    171       CALL wrk_alloc( jpi, jpj, env_up, env_dw, env0, env1, env2, env3 ) 
    172       CALL wrk_alloc( jpk, z_gsigw1, z_gsigt1, z_esigw1, z_esigt1 ) 
    173       ! 
    174       ! Initializing to 0.0 some arrays: 
    175       ! 
    176       ! *) gdept_1d, gdepw_1d : depth of T- and W-point (m) 
    177       ! *) e3t_1d  , e3w_1d   : scale factors at T- and W-levels (m) 
    178       ! *) z_gsigt1, z_gsigw1 : MES-coordinates (0.<= s <= 1.) of T- and W-points 
    179  
    180       gdept_1d = 0._wp   ;   gdepw_1d = 0._wp   ; 
    181       e3t_1d   = 0._wp   ;   e3w_1d   = 0._wp   ; 
    182       z_gsigw1 = 0._wp   ;   z_gsigt1 = 0._wp   ; 
    183       z_esigw1 = 0._wp   ;   z_esigt1 = 0._wp   ; 
    184  
    185158      ! Initialize mbathy to the maximum ocean level available 
    186  
    187159      mbathy(:,:) = jpkm1 
    188160 
    189161      ! Defining scosrf and scobot 
    190  
    191162      scosrf(:,:) = 0._wp             ! ocean surface depth (here zero: no under ice-shelf sea) 
    192163      scobot(:,:) = bathy(:,:)        ! ocean bottom  depth 
    193164 
    194       !================================================================== 
    195       ! 1. Defining 1D MEs-levels and computing their depths  
    196       !    (they will be used in diawri.F90) 
    197       !================================================================== 
    198  
    199       ! Computing gdept_1d and gdepw_1d 
    200  
    201       DO je = 1, tot_env ! LOOP over all the requested envelopes 
    202          IF (je == 1) THEN 
    203             s_1st  = 1 
    204             num_s  = nn_slev(je) 
    205             dep_up = 0._wp      ! surface 
    206          ELSE 
    207             s_1st = s_1st + num_s - 1 
    208             num_s  = nn_slev(je)  + 1 
    209             dep_up = rn_ebot_max(je-1) 
    210          ENDIF 
    211  
    212          dep_dw  = rn_ebot_max(je) 
    213          env_hc  = rn_e_hc(je) 
    214          coefa   = rn_e_th(je) 
    215          coefb   = rn_e_bb(je) 
    216          alpha   = rn_e_al(je) 
    217  
    218          IF (nn_strt(je) == 2) THEN 
    219             coefa = coefa / (dep_dw - dep_up) 
    220             coefb = coefb + ((dep_dw - dep_up)*rn_e_ba(je)) 
    221             coefb = 1.0_wp - (coefb/(dep_dw - dep_up)) 
    222          ENDIF 
    223  
    224          IF ( isodd(je) ) THEN 
    225             DO jk = 1, num_s 
    226                ind = s_1st+jk-1 
    227                ! W-GRID 
    228                z_gsigw1(ind) = -s_coord(sigma(1,jk,num_s), coefa, & 
    229                                         coefb, alpha, num_s, nn_strt(je)) 
    230                ! T-GRID 
    231                z_gsigt1(ind) = -s_coord(sigma(0,jk,num_s), coefa, & 
    232                                         coefb, alpha, num_s, nn_strt(je)) 
    233                zcoefw = z_gsigw1(ind) 
    234                zcoeft = z_gsigt1(ind) 
    235                !                   
    236                IF (nn_strt(je) /= 2) THEN 
    237                    gdept_1d(ind) = z_mes(dep_up, dep_dw, zcoeft, & 
    238                                         -sigma(0,jk,num_s), env_hc) 
    239                    gdepw_1d(ind) = z_mes(dep_up, dep_dw, zcoefw, & 
    240                                         -sigma(1,jk,num_s), env_hc) 
    241                ELSE 
    242                    gdept_1d(ind) = z_mes(dep_up, dep_dw, zcoeft, & 
    243                                         -sigma(0,jk,num_s), 0._wp) 
    244                    gdepw_1d(ind) = z_mes(dep_up, dep_dw, zcoefw, & 
    245                                         -sigma(1,jk,num_s), 0._wp) 
    246                END IF 
    247             ENDDO 
    248             ! 
    249          ELSE 
    250             ! Ensure that for splines we do not use SF12 
    251             !IF (nn_strt(je) > 1) nn_strt(je) = 1 
    252             ! 
    253             ! Interval's endpoints TAU are 0 and 1. 
    254             tau(1) = 0._wp 
    255             tau(n) = 1._wp 
    256  
    257             IF ( je == 2 ) THEN 
    258                dep0 = 0._wp 
    259             ELSE  
    260                dep0 = rn_ebot_max(je-2) 
    261             END IF  
    262             dep1   = dep_up 
    263             dep2   = dep_dw 
    264             IF ( je < tot_env ) dep3 = rn_ebot_max(je+1) 
    265             ! 
    266             ! Computing derivative analytically 
    267             ! 
    268             ! 1. Derivative for upper envelope 
    269             ibcbeg = 1 
    270             IF (nn_strt(je-1) == 2) THEN 
    271                alpha    = rn_e_al(je-1) 
    272                num_s_up = nn_slev(je-1) 
    273                coefa    = rn_e_th(je-1) / (dep1-dep0) 
    274                coefb    = rn_e_bb(je-1) + ((dep1-dep0)*rn_e_ba(je-1)) 
    275                coefb    = 1.0_wp-(coefb/(dep1-dep0)) 
    276             ELSE 
    277                alpha    = 0._wp 
    278                num_s_up = 1 
    279                coefa    = rn_e_th(je-1) 
    280                coefb    = rn_e_bb(je-1) 
    281             END IF 
    282             D1s_up = D1s_coord(-1._wp, coefa, coefb, & 
    283                                alpha, num_s_up, nn_strt(je-1)) 
    284             d(1,1) = dep_up 
    285             IF (nn_strt(je-1) == 2) THEN 
    286                d(2,1) = D1z_mes(dep0, dep1, D1s_up, 0._wp, 1) 
    287             ELSE 
    288                d(2,1) = D1z_mes(dep0, dep1, D1s_up, rn_e_hc(je-1), 1) 
    289             END IF 
    290             ! 
    291             ! 2. Derivative for deeper envelope 
    292             IF ( je < tot_env ) THEN 
    293                ibcend = 1 ! Bound. cond for the 1st derivative           
    294                IF (nn_strt(je+1) == 2) THEN 
    295                   alpha    = rn_e_al(je+1) 
    296                   num_s_dw = nn_slev(je+1) 
    297                   coefa    = rn_e_th(je+1) / (dep3-dep2) 
    298                   coefb    = rn_e_bb(je+1) + ((dep3-dep2)*rn_e_ba(je+1)) 
    299                   coefb    = 1.0_wp-(coefb/(dep3-dep2)) 
    300                ELSE 
    301                   alpha    = 0._wp 
    302                   num_s_dw = 1 
    303                   coefa    = rn_e_th(je+1) 
    304                   coefb    = rn_e_bb(je+1) 
    305                END IF 
    306                D1s_dw = D1s_coord(0._wp, coefa, coefb, & 
    307                                   alpha, num_s_dw, nn_strt(je+1)) 
    308                d(1,2) = dep_dw 
    309                IF (nn_strt(je+1) == 2) THEN 
    310                   d(2,2) = D1z_mes(dep2, dep3, D1s_dw, 0._wp, 1) 
    311                ELSE 
    312                   d(2,2) = D1z_mes(dep2, dep3, D1s_dw, rn_e_hc(je+1), 1) 
    313                END IF 
    314             ELSE 
    315                ibcend = 2     ! Bound. cond for the 2nd derivative 
    316                d(2,2) = 0._wp ! 2nd der = 0 
    317             ENDIF 
    318             ! 
    319             ! Computing spline's coefficients 
    320             IF ( lwp ) THEN 
    321                WRITE(numout,*) "BOUNDARY COND. to CONSTRAIN SPLINES:" 
    322                WRITE(numout,*) "ibcbeg: ", ibcbeg 
    323                WRITE(numout,*) "ibcend: ", ibcend 
    324             ENDIF 
    325             c = cub_spl(tau, d, n, ibcbeg, ibcend ) 
    326  
    327             IF ( lwp ) THEN 
    328                WRITE(numout,*) "Subzone between depth ", dep1, " and depth ", dep2  
    329                WRITE(numout,*) "" 
    330                WRITE(numout,*) "D1s_up = ", D1s_up 
    331                WRITE(numout,*) "1st deriv. of z(mes) at depth ", dep1, " is ", d(2,1) 
    332                WRITE(numout,*) "D1s_dw = ", D1s_dw 
    333                WRITE(numout,*) "1st deriv. of z(mes) at depth ", dep2, " is ", d(2,2) 
    334                WRITE(numout,*) "" 
    335                WRITE(numout,*) "CUBIC SPLINE COEFFICENTS" 
    336                WRITE(numout,*) "c1 = ", c(1,1) 
    337                WRITE(numout,*) "c2 = ", c(2,1) 
    338                WRITE(numout,*) "c3 = ", c(3,1) 
    339                WRITE(numout,*) "c4 = ", c(4,1) 
    340                WRITE(numout,*) "" 
    341                WRITE(numout,*) "CUB. SPL. at s = ", -sigma(1, 1, num_s),      & 
    342                                " is ", z_cub_spl(-sigma(1, 1, num_s), tau, c) 
    343                WRITE(numout,*) "1st deriv. of CUB. SPL. at s = ", -sigma(1, 1, num_s) ,     &  
    344                                " is ", D1z_cub_spl(-sigma(1, 1, num_s), tau, c, 1) 
    345                WRITE(numout,*) "CUB. SPL. at s = ", -sigma(1, num_s, num_s) , & 
    346                                " is ", z_cub_spl(-sigma(1, num_s, num_s), tau, c) 
    347                WRITE(numout,*) "1st deriv. of CUB. SPL. at s = ", -sigma(1, num_s, num_s) , & 
    348                                " is ", D1z_cub_spl(-sigma(1, num_s, num_s), tau, c, 1) 
    349             END IF 
    350  
    351             env_hc  = rn_e_hc(je) 
    352             coefa   = rn_e_th(je) 
    353             coefb   = rn_e_bb(je) 
    354             alpha   = rn_e_al(je) 
    355  
    356             DO jk = 1, num_s 
    357                ind = s_1st+jk-1 
    358                ! W-GRID 
    359                z_gsigw1(s_1st+jk-1) = -s_coord( sigma(1,jk,num_s), coefa, & 
    360                                                 coefb, alpha, num_s, nn_strt(je) ) 
    361                ! T-GRID 
    362                z_gsigt1(s_1st+jk-1) = -s_coord( sigma(0,jk,num_s), coefa, & 
    363                                                 coefb, alpha, num_s, nn_strt(je) ) 
    364                st  = z_gsigt1(ind) 
    365                sw  = z_gsigw1(ind) 
    366                gdept_1d(ind) = z_cub_spl(st, tau, c) 
    367                gdepw_1d(ind) = z_cub_spl(sw, tau, c) 
    368             END DO 
    369   
    370          END IF 
    371  
    372       END DO 
    373  
    374       !================================================================== 
    375       ! 2. Computing derivative of MEs-coordinate and scale factors 
    376       !    as finite differences 
    377       !================================================================== 
    378  
    379       ! Computing derivative of MEs-coordinate as finite differences 
    380       DO jk = 1, jpkm1 
    381          z_esigt1(jk)   = z_gsigw1(jk+1) - z_gsigw1(jk) 
    382          z_esigw1(jk+1) = z_gsigt1(jk+1) - z_gsigt1(jk) 
    383       END DO 
    384       ! Surface 
    385       z_esigw1( 1 ) = 2._wp * ( z_gsigt1( 1 ) - z_gsigw1( 1 ) ) 
    386       ! Bottom 
    387       z_esigt1(jpk) = 2._wp * ( z_gsigt1(jpk) - z_gsigw1(jpk) ) 
    388  
    389       ! Computing scale factors as finite differences 
    390       DO jk=1,jpkm1 
    391          e3t_1d(jk)   = gdepw_1d(jk+1) - gdepw_1d(jk) 
    392          e3w_1d(jk+1) = gdept_1d(jk+1) - gdept_1d(jk) 
    393       ENDDO 
    394       ! Bottom: 
    395       jk = jpk 
    396       e3t_1d(jk) = 2.0_wp * (gdept_1d(jk) - gdepw_1d(jk)) 
    397       ! and Surface  (H. Liu) 
    398       jk = 1 
    399       e3w_1d(jk) = 2.0_wp * (gdept_1d(1) - gdepw_1d(1)) 
    400   
    401       ! Control print 
    402       IF ( lwp ) THEN 
    403         WRITE(numout,*) "" 
    404         WRITE(numout,*) "1a. MEs-levels and first derivative" 
    405         WRITE(numout,*) "" 
    406         WRITE(numout,*) "      k     z_gsigw1      z_esigw1     z_gsigt1      z_esigt1" 
    407         WRITE(numout,*) "" 
    408         DO jk = 1, jpk 
    409            WRITE(numout,*) jk, z_gsigw1(jk), z_esigw1(jk), z_gsigt1(jk), z_esigt1(jk) 
    410         ENDDO 
    411         WRITE(numout,*) "-----------------------------------------------------------------" 
    412          
    413         WRITE(numout,*) "" 
    414         WRITE(numout,*) "2a. Depth of MEs-levels and scale factors" 
    415         WRITE(numout,*) "" 
    416         WRITE(numout,*) "      k     gdepw_1d      e3w_1d       gdepw_1d      e3w_1d" 
    417         WRITE(numout,*) "" 
    418         DO jk = 1, jpk 
    419            WRITE(numout,*) jk, gdepw_1d(jk), e3w_1d(jk), gdept_1d(jk), e3t_1d(jk) 
    420         ENDDO 
    421         WRITE(numout,*) "-----------------------------------------------------------------" 
    422       ENDIF 
    423  
    424       ! Checking monotonicity for cubic splines 
    425       DO jk = 1, jpk-1 
    426          IF ( gdept_1d(jk+1) < gdept_1d(jk) ) THEN 
    427             WRITE(ctlmes,*) 'NOT MONOTONIC gdept_1d: change envelopes' 
    428             CALL ctl_stop( ctlmes ) 
    429          END IF 
    430          IF ( gdepw_1d(jk+1) < gdepw_1d(jk) ) THEN 
    431             WRITE(ctlmes,*) 'NOT MONOTONIC gdepw_1d: change envelopes' 
    432             CALL ctl_stop( ctlmes ) 
    433          END IF 
    434       END DO 
    435  
    436165      !======================== 
    437       ! 3. Compute 3D MEs mesh 
     166      ! Compute 3D MEs mesh 
    438167      !======================== 
    439168 
     
    449178            env_up = envlt(:,:,je-1) 
    450179         ENDIF 
    451  
    452180         env_dw  = envlt(:,:,je) 
    453          env_hc  = rn_e_hc(je) 
    454          alpha   = rn_e_al(je) 
    455181 
    456182         IF ( isodd(je) ) THEN 
     183            ! 
     184            env_hc  = rn_e_hc(je) 
     185            coefa   = rn_e_th(je) 
     186            coefb   = rn_e_bb(je) 
     187            alpha   = rn_e_al(je) 
    457188            ! 
    458189            DO jk = 1, num_s 
    459190               ind = s_1st+jk-1 
     191               ! 
    460192               DO jj = 1,jpj 
    461193                  DO ji = 1,jpi 
    462                      IF (env_dw(ji,jj) < env_hc) THEN    
    463                         ! shallow water, uniform sigma 
    464                         zcoefw = -sigma(1, jk, num_s)  
     194                     ! 
     195                     ! ------------------------------------------- 
     196                     ! Computing MEs-coordinates (-1 <= MEs <= 0) 
     197                     ! 
     198                     ! shallow water, uniform sigma 
     199                     IF (env_dw(ji,jj) < env_hc) THEN 
     200                        ! W-GRID    
     201                        zcoefw = -sigma(1, jk, num_s) 
     202                        ! T-GRID  
    465203                        zcoeft = -sigma(0, jk, num_s) 
    466204                        !IF (nn_strt(je) == 2) THEN 
     
    471209                        !   zcoeft = zcoeft*(env_hc/(env_dw(ji,jj)-env_up(ji,jj))) 
    472210                        !ENDIF 
    473                      ELSE                               ! deep water, stretched s 
     211                     ! deep water, stretched s 
     212                     ELSE 
    474213                        IF (nn_strt(je) == 2) THEN 
    475214                           coefa = rn_e_th(je) /  (env_dw(ji,jj)-env_up(ji,jj)) 
    476215                           coefb = rn_e_bb(je) + ((env_dw(ji,jj)-env_up(ji,jj))*rn_e_ba(je)) 
    477216                           coefb = 1.0_wp-(coefb/(env_dw(ji,jj)-env_up(ji,jj))) 
    478                            ! W-GRID 
    479                            zcoefw = -s_coord( sigma(1,jk,num_s), coefa, & 
    480                                               coefb, alpha, num_s, nn_strt(je) ) 
    481                            ! T-GRID 
    482                            zcoeft = -s_coord( sigma(0,jk,num_s), coefa, & 
    483                                               coefb, alpha, num_s, nn_strt(je) ) 
    484                         ELSE 
    485                            zcoefw = z_gsigw1(ind) 
    486                            zcoeft = z_gsigt1(ind) 
    487217                        END IF 
     218                        ! W-GRID 
     219                        zcoefw = -s_coord( sigma(1,jk,num_s), coefa, & 
     220                                           coefb, alpha, num_s, nn_strt(je) ) 
     221                        ! T-GRID 
     222                        zcoeft = -s_coord( sigma(0,jk,num_s), coefa, & 
     223                                           coefb, alpha, num_s, nn_strt(je) ) 
    488224                     ENDIF 
    489225                     ! 
     226                     z_gsigw3(ji,jj,ind) = zcoefw 
     227                     z_gsigt3(ji,jj,ind) = zcoeft 
     228                     ! 
     229                     ! -------------------------------------------------- 
     230                     ! Computing model levels depths 
    490231                     IF (nn_strt(je) /= 2) THEN 
    491232                        gdept_0(ji,jj,ind) = z_mes(env_up(ji,jj), env_dw(ji,jj), zcoeft, & 
     
    591332                  !   WRITE(numout,*) "ibcend: ", ibcend 
    592333                  !ENDIF 
     334 
    593335                  c = cub_spl(tau, d, n, ibcbeg, ibcend ) 
    594336 
    595                   ! Computing model level depths 
     337                  env_hc  = rn_e_hc(je) 
     338                  coefa   = rn_e_th(je) 
     339                  coefb   = rn_e_bb(je) 
     340                  alpha   = rn_e_al(je) 
     341 
    596342                  DO jk = 1, num_s 
    597343                     ind = s_1st+jk-1 
    598                      st  = z_gsigt1(ind) 
    599                      sw  = z_gsigw1(ind) 
    600                      gdept_0(ji,jj,ind) = z_cub_spl(st, tau, c) 
    601                      gdepw_0(ji,jj,ind) = z_cub_spl(sw, tau, c) 
     344                     ! 
     345                     ! Computing stretched distribution of interpolation points 
     346                     ! W-GRID 
     347                     zcoefw = -s_coord( sigma(1,jk,num_s), coefa, & 
     348                                        coefb, alpha, num_s, nn_strt(je) ) 
     349                     ! T-GRID 
     350                     zcoeft = -s_coord( sigma(0,jk,num_s), coefa, & 
     351                                        coefb, alpha, num_s, nn_strt(je) )                      
     352                     ! 
     353                     z_gsigw3(ji,jj,ind) = zcoefw 
     354                     z_gsigt3(ji,jj,ind) = zcoeft 
     355                     ! 
     356                     gdept_0(ji,jj,ind) = z_cub_spl(zcoeft, tau, c) 
     357                     gdepw_0(ji,jj,ind) = z_cub_spl(zcoefw, tau, c) 
    602358                  END DO 
    603359 
    604                   ! Control print 
    605                   IF ( lwp ) THEN 
    606                      IF ( jj == 25 .AND. ji == 7 ) THEN 
    607                         WRITE(numout,*) "" 
    608                   !   WRITE(numout,*) "env_up = ", env1(ji,jj) 
    609                         WRITE(numout,*) "1st deriv. at TAU = 0._wp is ", d(2,1) 
    610                   !   WRITE(numout,*) "env_dw = ", env2(ji,jj) 
    611                         WRITE(numout,*) "1st deriv. at TAU = 1._wp is ", d(2,2) 
    612                   !   WRITE(numout,*) "CUBIC SPLINE at s = ", -sigma(1, 1, num_s) , " is ", & 
    613                   !                   z_cub_spl(-sigma(1, 1, num_s), tau, c) 
    614                         WRITE(numout,*) "1st deriv. at s = ", -sigma(1, 1, num_s) , " is ", & 
    615                                         D1z_cub_spl(-sigma(1, 1, num_s), tau, c, 1) 
    616                   !   WRITE(numout,*) "CUBIC SPLINE at s = ", -sigma(1, num_s, num_s) , " is ", & 
    617                   !                   z_cub_spl(-sigma(1, num_s, num_s), tau, c) 
    618                         WRITE(numout,*) "1st deriv. at s = ", -sigma(1, num_s, num_s) , " is ", & 
    619                                         D1z_cub_spl(-sigma(1, num_s, num_s), tau, c, 1) 
    620                      END IF 
    621                   END IF 
    622   
    623360               END DO 
    624361            END DO 
     
    644381      jk = jpk 
    645382      e3t_0(:,:,jk) = 2.0_wp * (gdept_0(:,:,jk) - gdepw_0(:,:,jk)) 
    646  
    647       IF ( lwp ) THEN 
    648          jj = 25  
    649          ji = 7 
    650          WRITE(numout,*) "" 
    651          WRITE(numout,*) "      k    gdepw_0(k)    gdept_0(k)    e3w_0(k)    e3t_0(k)" 
    652          WRITE(numout,*) "" 
    653          DO jk = 1, jpk 
    654             WRITE(numout,*) jk, gdepw_0(ji,jj,jk), gdept_0(ji,jj,jk), e3w_0(ji,jj,jk), e3t_0(ji,jj,jk) 
    655          ENDDO 
    656       END IF 
    657  
    658       !============================== 
    659       ! Computing gdept3w 
    660       !============================== 
    661  
    662       gde3w_0(:,:,1) = 0.5 * e3w_0(:,:,1) 
    663       DO jk = 2, jpk 
    664          gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk) 
    665       END DO 
    666       e3t_0(:,:,jpk) = 2.0_wp * (gdept_0(:,:,jpk) - gdepw_0(:,:,jpk)) 
    667383 
    668384      !============================== 
     
    825541      !================================================================================ 
    826542 
     543      ! Extracting MEs depth profile in the shallowest point of the deepest  
     544      ! envelope for a first check of monotonicity of transformation  
     545      ! (also useful for debugging) 
     546      pmsk(:,:) = 0.0 
     547      WHERE ( bathy > 0.0 ) pmsk = 1.0 
     548      CALL mpp_minloc( envlt(:,:,tot_env), pmsk(:,:), mindep, eii, eij ) 
     549      IF ((mi0(eii)>1 .AND. mi0(eii)<jpi) .AND. (mj0(eij)>1 .AND. mj0(eij)<jpj)) THEN 
     550         irnk = mpprank 
     551         DO je = 1, tot_env 
     552            rn_ebot_max(je) = envlt(mi0(eii),mj0(eij),je) 
     553         END DO 
     554         z_gsigt1(:) = z_gsigt3(mi0(eii),mj0(eij),:) 
     555         z_gsigw1(:) = z_gsigw3(mi0(eii),mj0(eij),:) 
     556         gdept1(:)   = gdept_0(mi0(eii),mj0(eij),:) 
     557         gdepw1(:)   = gdepw_0(mi0(eii),mj0(eij),:) 
     558         e3t1(:)     = e3t_0(mi0(eii),mj0(eij),:) 
     559         e3w1(:)     = e3w_0(mi0(eii),mj0(eij),:) 
     560      ELSE 
     561         irnk = -1 
     562      END IF 
     563      IF( lk_mpp ) CALL mppsync 
     564      IF( lk_mpp ) CALL mpp_max(irnk) 
     565      IF( lk_mpp ) CALL mppbcast_a_real(rn_ebot_max, max_nn_env, irnk ) 
     566      IF( lk_mpp ) CALL mppbcast_a_real(z_gsigt1   , jpk       , irnk ) 
     567      IF( lk_mpp ) CALL mppbcast_a_real(z_gsigw1   , jpk       , irnk ) 
     568      IF( lk_mpp ) CALL mppbcast_a_real(gdept1     , jpk       , irnk ) 
     569      IF( lk_mpp ) CALL mppbcast_a_real(gdepw1     , jpk       , irnk ) 
     570      IF( lk_mpp ) CALL mppbcast_a_real(e3t1       , jpk       , irnk ) 
     571      IF( lk_mpp ) CALL mppbcast_a_real(e3w1       , jpk       , irnk ) 
     572 
     573      IF( lwp ) THEN 
     574        WRITE(numout,*) "" 
     575        WRITE(numout,*) "mes_build:" 
     576        WRITE(numout,*) "~~~~~~~~~" 
     577        WRITE(numout,*) "" 
     578        WRITE(numout,*) " FIRST CHECK: Checking MEs-levels profile in the shallowest point of the last envelope:" 
     579        WRITE(numout,*) "              it is the most likely point where monotonicty of splines may be violeted. " 
     580        WRITE(numout,*) "" 
     581        DO je = 1, tot_env 
     582           WRITE(numout,*) '              * depth of envelope ', je, ' at point (',eii,',',eij,') is ', rn_ebot_max(je) 
     583        END DO 
     584        WRITE(numout,*) "" 
     585        WRITE(numout,*) "      MEs-coordinates" 
     586        WRITE(numout,*) "" 
     587        WRITE(numout,*) "      k     z_gsigw1      z_gsigt1" 
     588        WRITE(numout,*) "" 
     589        DO jk = 1, jpk 
     590           WRITE(numout,*) '   ', jk, z_gsigw1(jk), z_gsigt1(jk) 
     591        ENDDO 
     592        WRITE(numout,*) "" 
     593        WRITE(numout,*) "-----------------------------------------------------------------" 
     594        WRITE(numout,*) "" 
     595        WRITE(numout,*) "      MEs-levels depths and scale factors" 
     596        WRITE(numout,*) "" 
     597        WRITE(numout,*) "      k     gdepw1      e3w1       gdept1      e3t1" 
     598        WRITE(numout,*) "" 
     599        DO jk = 1, jpk  
     600           WRITE(numout,*) '   ', jk, gdepw1(jk), e3w1(jk), gdept1(jk), e3t1(jk) 
     601        ENDDO 
     602        WRITE(numout,*) "-----------------------------------------------------------------" 
     603      ENDIF 
     604 
     605      ! Checking monotonicity for cubic splines 
     606      DO jk = 1, jpk-1 
     607         IF ( gdept1(jk+1) < gdept1(jk) ) THEN 
     608            WRITE(ctlmes,*) 'NOT MONOTONIC gdept_0: change envelopes' 
     609            CALL ctl_stop( ctlmes ) 
     610         END IF 
     611         IF ( gdepw1(jk+1) < gdepw1(jk) ) THEN 
     612            WRITE(ctlmes,*) 'NOT MONOTONIC gdepw_0: change envelopes' 
     613            CALL ctl_stop( ctlmes ) 
     614         END IF 
     615         IF ( e3t1(jk) < 0.0 ) THEN 
     616            WRITE(ctlmes,*) 'NEGATIVE e3t_0: change envelopes' 
     617            CALL ctl_stop( ctlmes ) 
     618         END IF 
     619         IF ( e3w1(jk) < 0.0 ) THEN 
     620            WRITE(ctlmes,*) 'NEGATIVE e3w_0: change envelopes' 
     621            CALL ctl_stop( ctlmes ) 
     622         END IF 
     623      END DO 
     624 
     625      ! CHECKING THE WHOLE DOMAIN 
    827626      DO ji = 1, jpi 
    828627         DO jj = 1, jpj 
    829             DO jk = 1, mbathy(ji,jj) 
     628            DO jk = 1, jpk !mbathy(ji,jj) 
    830629               ! check coordinate is monotonically increasing 
    831630               IF (e3w_0(ji,jj,jk) <= 0._wp .OR. e3t_0(ji,jj,jk) <= 0._wp ) THEN 
    832                   WRITE(ctmp1,*) 'ERROR zgr_mes :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
    833                   WRITE(numout,*) 'ERROR zgr_mes :   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
     631                  WRITE(ctmp1,*) 'ERROR mes_build:   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
     632                  WRITE(numout,*) 'ERROR mes_build:   e3w   or e3t   =< 0  at point (i,j,k)= ', ji, jj, jk 
    834633                  WRITE(numout,*) 'e3w',e3w_0(ji,jj,:) 
    835634                  WRITE(numout,*) 'e3t',e3t_0(ji,jj,:) 
     
    838637               ! and check it has never gone negative 
    839638               IF ( gdepw_0(ji,jj,jk) < 0._wp .OR. gdept_0(ji,jj,jk) < 0._wp ) THEN 
    840                   WRITE(ctmp1,*) 'ERROR zgr_mes :   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk 
    841                   WRITE(numout,*) 'ERROR zgr_mes :   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk 
     639                  WRITE(ctmp1,*) 'ERROR mes_build:   gdepw or gdept =< 0  at point (i,j,k)= ', ji, jj, jk 
     640                  WRITE(numout,*) 'ERROR mes_build:   gdepw   or gdept   =< 0  at point (i,j,k)= ', ji, jj, jk 
    842641                  WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:) 
    843642                  WRITE(numout,*) 'gdept',gdept_0(ji,jj,:) 
    844643                  CALL ctl_stop( ctmp1 ) 
    845644               ENDIF 
    846 !              ! and check it never exceeds the total depth 
    847 !               IF ( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN 
    848 !                  WRITE(ctmp1,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk 
    849 !                  WRITE(numout,*) 'ERROR zgr_sco :   gdepw > hbatt  at point (i,j,k)= ', ji, jj, jk 
    850 !                  WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:) 
    851 !                  CALL ctl_stop( ctmp1 ) 
    852 !               ENDIF 
    853645            END DO 
    854  
    855 !            DO jk = 1, mbathy(ji,jj)-1 
    856 !               ! and check it never exceeds the total depth 
    857 !               IF ( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN 
    858 !                  WRITE(ctmp1,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk 
    859 !                  WRITE(numout,*) 'ERROR zgr_sco :   gdept > hbatt  at point (i,j,k)= ', ji, jj, jk 
    860 !                  WRITE(numout,*) 'gdept',fsdept(ji,jj,:) 
    861 !                  CALL ctl_stop( ctmp1 ) 
    862 !               ENDIF 
    863 !            END DO 
    864646         END DO 
    865647      END DO 
    866648      ! 
    867649      CALL wrk_dealloc( jpi, jpj, max_nn_env , envlt ) 
    868       CALL wrk_dealloc( jpi, jpj, env_up, env_dw ) 
    869       CALL wrk_dealloc( jpk, z_gsigw1, z_gsigt1, z_esigw1, z_esigt1 ) 
    870650      ! 
    871651      IF( nn_timing == 1 )  CALL timing_stop('mes_build') 
     
    876656 
    877657   SUBROUTINE mes_init 
     658      !!----------------------------------------------------------------------------- 
     659      !!                  ***  ROUTINE mes_init  *** 
     660      !!                      
     661      !! ** Purpose : Initilaise a Multi Enveloped s-coordinate (MEs) system 
     662      !! 
     663      !!----------------------------------------------------------------------------- 
     664 
    878665      CHARACTER(lc)                :: env_name   ! name of the externally defined envelope 
    879       INTEGER                      :: ji, jj, je, eii, eij 
    880       INTEGER                      :: cor_env, ios, inum 
    881       REAL(wp)                     :: rn_bot_min ! minimum depth of ocean bottom (>0) (m) 
    882       REAL(wp)                     :: rn_bot_max ! maximum depth of ocean bottom (= ocean depth) (>0) (m) 
    883       REAL(wp)                     :: mindep 
    884       REAL(wp), DIMENSION(jpi,jpj) :: pmsk, ewrk ! working arrays 
     666      INTEGER                      :: ji, jj, je, inum 
     667      INTEGER                      :: cor_env, ios 
     668      REAL(wp)                     :: rn_bot_min ! min depth of ocean bottom (>0) (m) 
     669      REAL(wp)                     :: rn_bot_max ! max depth of ocean bottom (= ocean depth) (>0) (m) 
    885670 
    886671      NAMELIST/namzgr_mes/rn_bot_min , rn_bot_max , ln_envl, & 
     
    893678      ! 
    894679      CALL wrk_alloc( jpi, jpj, max_nn_env , envlt ) 
    895       pmsk(:,:) = 0.0 
    896680      ! 
    897681      ! Initialising some arrays 
    898       rn_ebot_max(:) = 0.0 
    899682      max_dep_env(:) = 0.0 
    900683      min_dep_env(:) = 0.0 
     
    963746         IF (MAXVAL(envlt(:,:,je+1)) < MAXVAL(envlt(:,:,je))) CALL ctl_stop( ctlmes ) 
    964747      ENDDO 
    965       ! 5) Computing max depth of envelopes in the shallowest 
    966       !    point of the deeper envelope for first check of monotonicity 
    967       !    of transformation and computing 1D MEs-levels depths 
    968       !    (used by diawri.F90) 
    969       WHERE ( bathy > 0.0 ) pmsk = 1.0 
    970       ewrk(:,:) = envlt(:,:,tot_env) 
    971       CALL mpp_minloc( ewrk(:,:), pmsk(:,:), mindep, eii, eij ) 
    972       IF ((mi0(eii)>1 .AND. mi0(eii)<jpi) .AND. (mj0(eij)>1 .AND. mj0(eij)<jpj)) THEN         
    973          inum = mpprank 
    974          DO je = 1, tot_env 
    975             rn_ebot_max(je) = envlt(mi0(eii),mj0(eij),je) 
    976          END DO 
    977       ELSE 
    978          inum = -1 
    979       END IF 
    980       IF( lk_mpp ) CALL mppsync 
    981       IF( lk_mpp ) CALL mpp_max(inum) 
    982       IF( lk_mpp ) CALL mppbcast_a_real(rn_ebot_max, max_nn_env, inum ) 
     748      ! 5) Computing max and min depths of envelopes 
    983749      DO je = 1, tot_env 
    984750         max_dep_env(je) = MAXVAL(envlt(:,:,je)) 
Note: See TracChangeset for help on using the changeset viewer.