[15126] | 1 | MODULE mes |
---|
| 2 | !!============================================================================== |
---|
[15271] | 3 | !! *** MODULE mes *** |
---|
[15126] | 4 | !! Ocean initialization : Multiple Enveloped s coordinate (MES) |
---|
| 5 | !!============================================================================== |
---|
[15134] | 6 | !! NEMO 3.6 ! 2017-05 D. Bruciaferri |
---|
| 7 | !! NEMO 4.0.4 ! 2021-07 D. Bruciaferri |
---|
[15126] | 8 | !!---------------------------------------------------------------------- |
---|
[15134] | 9 | ! |
---|
[15126] | 10 | USE oce ! ocean variables |
---|
| 11 | USE dom_oce ! ocean domain |
---|
| 12 | USE closea ! closed seas |
---|
| 13 | ! |
---|
| 14 | USE in_out_manager ! I/O manager |
---|
| 15 | USE iom ! I/O library |
---|
| 16 | USE lbclnk ! ocean lateral boundary conditions (or mpp link) |
---|
| 17 | USE lib_mpp ! distributed memory computing library |
---|
| 18 | USE wrk_nemo ! Memory allocation |
---|
| 19 | USE timing ! Timing |
---|
| 20 | |
---|
| 21 | IMPLICIT NONE |
---|
| 22 | PRIVATE |
---|
| 23 | |
---|
[15129] | 24 | PUBLIC mes_build ! called by zgrmes.F90 |
---|
| 25 | ! |
---|
| 26 | ! Envelopes and stretching parameters |
---|
[15126] | 27 | CHARACTER(lc) :: ctlmes ! control message error (lc=256) |
---|
| 28 | INTEGER, PARAMETER :: max_nn_env = 5 ! Maximum allowed number of envelopes. |
---|
[15129] | 29 | INTEGER :: tot_env ! Tot number of requested envelopes |
---|
[15126] | 30 | LOGICAL :: ln_envl(max_nn_env) ! Array with flags (T/F) specifing which envelope is used |
---|
| 31 | INTEGER :: nn_strt(max_nn_env) ! Array specifing the stretching function for each envelope: |
---|
| 32 | ! Madec 1996 (0), Song and Haidvogel 1994 (1) |
---|
[15131] | 33 | REAL(wp) :: max_dep_env(max_nn_env) ! global maximum depth of envelopes |
---|
| 34 | REAL(wp) :: min_dep_env(max_nn_env) ! global minimum depth of envelopes |
---|
[15126] | 35 | INTEGER :: nn_slev(max_nn_env) ! Array specifing number of levels of each enveloped vertical zone |
---|
| 36 | REAL(wp) :: rn_e_hc(max_nn_env) ! Array specifing critical depth for transition to stretched |
---|
| 37 | ! coordinates of each envelope |
---|
| 38 | REAL(wp) :: rn_e_th(max_nn_env) ! Array specifing surface control parameter (0<=th<=20) of SH94 |
---|
| 39 | ! or rn_zs of SF12 for each vertical sub-zone |
---|
| 40 | REAL(wp) :: rn_e_bb(max_nn_env) ! Array specifing bottom control parameter (0<=bb<=1) of SH94 |
---|
| 41 | ! or rn_zb_b of SF12 for each vertical sub-zone |
---|
| 42 | REAL(wp) :: rn_e_ba(max_nn_env) ! Array specifing bottom control parameter rn_zb_a of SF12 for |
---|
| 43 | ! each vertical sub-zone |
---|
| 44 | REAL(wp) :: rn_e_al(max_nn_env) ! Array specifing stretching parameter rn_alpha of SF12 for |
---|
| 45 | ! each vertical sub-zone |
---|
[15278] | 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.) |
---|
[15129] | 51 | ! |
---|
| 52 | REAL(wp), POINTER, DIMENSION(:,:,:) :: envlt ! array for the envelopes |
---|
[15126] | 53 | |
---|
| 54 | !! * Substitutions |
---|
| 55 | # include "vectopt_loop_substitute.h90" |
---|
| 56 | |
---|
| 57 | CONTAINS |
---|
| 58 | |
---|
| 59 | ! ===================================================================================================== |
---|
| 60 | |
---|
[15129] | 61 | SUBROUTINE mes_build |
---|
[15126] | 62 | !!----------------------------------------------------------------------------- |
---|
[15129] | 63 | !! *** ROUTINE mes_build *** |
---|
[15126] | 64 | !! |
---|
| 65 | !! ** Purpose : define the Multi Enveloped S-coordinate (MES) system |
---|
| 66 | !! |
---|
| 67 | !! ** Method : s-coordinate defined respect to multiple |
---|
[15134] | 68 | !! externally defined envelope as detailed in |
---|
[15126] | 69 | !! |
---|
[15134] | 70 | !! Bruciaferri, Shapiro, Wobus, 2018. Oce. Dyn. |
---|
| 71 | !! https://doi.org/10.1007/s10236-018-1189-x |
---|
[15126] | 72 | !! |
---|
[15134] | 73 | !! Three options for stretching are given: |
---|
[15126] | 74 | !! |
---|
[15134] | 75 | !! *) nn_strt = 0 : Madec et al 1996 cosh/tanh function |
---|
| 76 | !! *) nn_strt = 1 : Song and Haidvogel 1994 sinh/tanh function |
---|
| 77 | !! *) nn_strt = 2 : Siddorn and Furner gamma function |
---|
[15126] | 78 | !! |
---|
| 79 | !!---------------------------------------------------------------------- |
---|
[15134] | 80 | !! SKETCH of the GEOMETRY OF A MEs-COORDINATE SYSTEM |
---|
[15126] | 81 | !! |
---|
| 82 | !! Lines represent W-levels: |
---|
| 83 | !! |
---|
| 84 | !! 0: First s-level part. The total number |
---|
| 85 | !! of levels in this part is controlled |
---|
| 86 | !! by the nn_slev(1) namelist parameter. |
---|
| 87 | !! |
---|
| 88 | !! Depth first W-lev: 0 m (surfcace) |
---|
| 89 | !! Depth last W-lev: depth of 1st envelope |
---|
| 90 | !! |
---|
| 91 | !! o: Second s-level part. The total number |
---|
| 92 | !! of levels in this part is controlled |
---|
| 93 | !! by the nn_slev(2) namelist parameter. |
---|
| 94 | !! |
---|
| 95 | !! Depth last W-lev: depth of 2nd envelope |
---|
| 96 | !! |
---|
| 97 | !! @: Third s-level part. The total number |
---|
| 98 | !! of levels in this part is controlled |
---|
| 99 | !! by the nn_slev(3) namelist parameter. |
---|
| 100 | !! |
---|
| 101 | !! Depth last W-lev: depth of 3rd envelope |
---|
| 102 | !! |
---|
| 103 | !! z |~~~~~~~~~~~~~~~~~~~~~~~0~~~~~~~~~~~~~~~~~~~~~~~ SURFACE |
---|
| 104 | !! | |
---|
| 105 | !! |-----------------------0----------------------- ln_envl(1) = true, nn_slev(1) = 3 |
---|
| 106 | !! | |
---|
| 107 | !! |=======================0======================= ENVELOPE 1 |
---|
| 108 | !! | |
---|
| 109 | !! |-----------------------o----------------------- |
---|
| 110 | !! | ln_envl(2) = true, nn_slev(2) = 3 |
---|
| 111 | !! |-----------------------o----------------------- |
---|
| 112 | !! | |
---|
| 113 | !! |¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬o¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬ ENVELOPE 2 |
---|
| 114 | !! | |
---|
| 115 | !! |-----------------------@----------------------- ln_envl(3) = true, nn_slev(3) = 2 |
---|
| 116 | !! | |
---|
| 117 | !! |_______________________@_______________________ ENVELOPE 3 |
---|
| 118 | !! \|/ |
---|
| 119 | !! |
---|
| 120 | !! |
---|
| 121 | !!----------------------------------------------------------------------------- |
---|
| 122 | !! Author: Diego Bruciaferri (University of Plymouth), September 2017 |
---|
| 123 | !!----------------------------------------------------------------------------- |
---|
| 124 | |
---|
| 125 | ! |
---|
[15134] | 126 | INTEGER :: ji, jj, jk, jl, je ! dummy loop argument |
---|
| 127 | INTEGER :: eii, eij, irnk ! dummy loop argument |
---|
| 128 | INTEGER :: num_s, s_1st, ind ! for loops over envelopes |
---|
| 129 | INTEGER :: num_s_up, num_s_dw ! for loops over envelopes |
---|
| 130 | REAL(wp) :: D1s_up, D1s_dw ! for loops over envelopes |
---|
| 131 | INTEGER :: ibcbeg, ibcend ! boundary condition for cubic splines |
---|
| 132 | INTEGER, PARAMETER :: n = 2 ! number of considered points for each |
---|
| 133 | ! ! interval |
---|
| 134 | REAL(wp) :: c(4,n) ! matrix for cubic splines coefficents |
---|
| 135 | REAL(wp) :: d(2,2) ! matrix for depth values and depth |
---|
| 136 | ! ! derivative at intervals' boundaries |
---|
| 137 | REAL(wp) :: tau(n) ! abscissas of intervals' boundaries |
---|
| 138 | REAL(wp) :: coefa, coefb ! for loops over envelopes |
---|
| 139 | REAL(wp) :: alpha, env_hc ! for loops over envelopes |
---|
| 140 | REAL(wp) :: zcoeft, zcoefw ! temporary scalars |
---|
| 141 | REAL(wp) :: max_env_up, min_env_up ! to identify geopotential levels |
---|
| 142 | REAL(wp) :: max_env_dw, min_env_dw ! to identify geopotential levels |
---|
| 143 | REAL(wp) :: mindep |
---|
| 144 | ! |
---|
| 145 | REAL(wp), DIMENSION(max_nn_env) :: rn_ebot_max |
---|
| 146 | REAL(wp), DIMENSION(jpk) :: z_gsigw1, z_gsigt1 ! 1D arrays for debugging |
---|
| 147 | REAL(wp), DIMENSION(jpk) :: gdepw1, gdept1 ! 1D arrays for debugging |
---|
| 148 | REAL(wp), DIMENSION(jpk) :: e3t1, e3w1 ! 1D arrays for debugging |
---|
| 149 | ! |
---|
| 150 | REAL(wp), DIMENSION(jpi,jpj) :: env_up, env_dw ! for loops over envelopes |
---|
| 151 | REAL(wp), DIMENSION(jpi,jpj) :: env0, env1, env2, env3 ! for loops over envelopes |
---|
| 152 | ! ! for cubic splines |
---|
| 153 | REAL(wp), DIMENSION(jpi,jpj) :: pmsk ! working array |
---|
| 154 | ! |
---|
| 155 | REAL(wp), DIMENSION(jpi,jpj,jpk) :: z_gsigw3, z_gsigt3 ! 3D arrays for debugging |
---|
[15126] | 156 | !!---------------------------------------------------------------------- |
---|
| 157 | ! |
---|
[15129] | 158 | IF( nn_timing == 1 ) CALL timing_start('mes_build') |
---|
[15126] | 159 | ! |
---|
[15129] | 160 | CALL mes_init |
---|
| 161 | ! |
---|
[15126] | 162 | ! Initialize mbathy to the maximum ocean level available |
---|
| 163 | mbathy(:,:) = jpkm1 |
---|
| 164 | |
---|
| 165 | ! Defining scosrf and scobot |
---|
| 166 | scosrf(:,:) = 0._wp ! ocean surface depth (here zero: no under ice-shelf sea) |
---|
| 167 | scobot(:,:) = bathy(:,:) ! ocean bottom depth |
---|
| 168 | |
---|
| 169 | !======================== |
---|
[15134] | 170 | ! Compute 3D MEs mesh |
---|
[15126] | 171 | !======================== |
---|
| 172 | |
---|
| 173 | DO je = 1, tot_env ! LOOP over all the requested envelopes |
---|
| 174 | |
---|
| 175 | IF (je == 1) THEN |
---|
| 176 | s_1st = 1 |
---|
| 177 | num_s = nn_slev(je) |
---|
| 178 | env_up = scosrf ! surface |
---|
| 179 | ELSE |
---|
| 180 | s_1st = s_1st + num_s - 1 |
---|
| 181 | num_s = nn_slev(je) + 1 |
---|
| 182 | env_up = envlt(:,:,je-1) |
---|
| 183 | ENDIF |
---|
| 184 | env_dw = envlt(:,:,je) |
---|
| 185 | |
---|
| 186 | IF ( isodd(je) ) THEN |
---|
[15131] | 187 | ! |
---|
[15134] | 188 | env_hc = rn_e_hc(je) |
---|
| 189 | coefa = rn_e_th(je) |
---|
| 190 | coefb = rn_e_bb(je) |
---|
| 191 | alpha = rn_e_al(je) |
---|
| 192 | ! |
---|
[15131] | 193 | DO jk = 1, num_s |
---|
| 194 | ind = s_1st+jk-1 |
---|
[15134] | 195 | ! |
---|
[15131] | 196 | DO jj = 1,jpj |
---|
| 197 | DO ji = 1,jpi |
---|
[15134] | 198 | ! |
---|
| 199 | ! ------------------------------------------- |
---|
| 200 | ! Computing MEs-coordinates (-1 <= MEs <= 0) |
---|
| 201 | ! |
---|
| 202 | ! shallow water, uniform sigma |
---|
| 203 | IF (env_dw(ji,jj) < env_hc) THEN |
---|
| 204 | ! W-GRID |
---|
| 205 | zcoefw = -sigma(1, jk, num_s) |
---|
| 206 | ! T-GRID |
---|
[15131] | 207 | zcoeft = -sigma(0, jk, num_s) |
---|
| 208 | !IF (nn_strt(je) == 2) THEN |
---|
| 209 | ! ! shallow water, z coordinates |
---|
| 210 | ! ! SF12 in AMM* configurations uses |
---|
| 211 | ! ! this (ln_sigcrit=.false.) |
---|
| 212 | ! zcoefw = zcoefw*(env_hc/(env_dw(ji,jj)-env_up(ji,jj))) |
---|
| 213 | ! zcoeft = zcoeft*(env_hc/(env_dw(ji,jj)-env_up(ji,jj))) |
---|
| 214 | !ENDIF |
---|
[15134] | 215 | ! deep water, stretched s |
---|
| 216 | ELSE |
---|
[15131] | 217 | IF (nn_strt(je) == 2) THEN |
---|
| 218 | coefa = rn_e_th(je) / (env_dw(ji,jj)-env_up(ji,jj)) |
---|
| 219 | coefb = rn_e_bb(je) + ((env_dw(ji,jj)-env_up(ji,jj))*rn_e_ba(je)) |
---|
| 220 | coefb = 1.0_wp-(coefb/(env_dw(ji,jj)-env_up(ji,jj))) |
---|
[15126] | 221 | END IF |
---|
[15134] | 222 | ! W-GRID |
---|
| 223 | zcoefw = -s_coord( sigma(1,jk,num_s), coefa, & |
---|
| 224 | coefb, alpha, num_s, nn_strt(je) ) |
---|
| 225 | ! T-GRID |
---|
| 226 | zcoeft = -s_coord( sigma(0,jk,num_s), coefa, & |
---|
| 227 | coefb, alpha, num_s, nn_strt(je) ) |
---|
[15131] | 228 | ENDIF |
---|
| 229 | ! |
---|
[15134] | 230 | z_gsigw3(ji,jj,ind) = zcoefw |
---|
| 231 | z_gsigt3(ji,jj,ind) = zcoeft |
---|
| 232 | ! |
---|
| 233 | ! -------------------------------------------------- |
---|
| 234 | ! Computing model levels depths |
---|
[15131] | 235 | IF (nn_strt(je) /= 2) THEN |
---|
| 236 | gdept_0(ji,jj,ind) = z_mes(env_up(ji,jj), env_dw(ji,jj), zcoeft, & |
---|
| 237 | -sigma(0, jk, num_s), env_hc) |
---|
| 238 | |
---|
| 239 | gdepw_0(ji,jj,ind) = z_mes(env_up(ji,jj), env_dw(ji,jj), zcoefw, & |
---|
| 240 | -sigma(1, jk, num_s), env_hc) |
---|
| 241 | ELSE |
---|
| 242 | gdept_0(ji,jj,ind) = z_mes(env_up(ji,jj), env_dw(ji,jj), zcoeft, & |
---|
| 243 | -sigma(0, jk, num_s), 0._wp) |
---|
| 244 | |
---|
| 245 | gdepw_0(ji,jj,ind) = z_mes(env_up(ji,jj), env_dw(ji,jj), zcoefw, & |
---|
| 246 | -sigma(1, jk, num_s), 0._wp) |
---|
| 247 | END IF |
---|
[15126] | 248 | END DO |
---|
| 249 | END DO |
---|
[15131] | 250 | END DO |
---|
| 251 | ! |
---|
[15126] | 252 | ELSE |
---|
[15131] | 253 | ! |
---|
[15126] | 254 | ! Interval's endpoints TAU are 0 and 1. |
---|
| 255 | tau(1) = 0._wp |
---|
| 256 | tau(n) = 1._wp |
---|
[15131] | 257 | ! |
---|
[15126] | 258 | IF ( je == 2 ) THEN |
---|
| 259 | env0 = scosrf |
---|
| 260 | ELSE |
---|
| 261 | env0 = envlt(:,:,je-2) |
---|
| 262 | END IF |
---|
| 263 | env1 = env_up |
---|
| 264 | env2 = env_dw |
---|
[15131] | 265 | IF ( je < tot_env ) env3 = envlt(:,:,je+1) |
---|
| 266 | ! |
---|
[15126] | 267 | DO jj = 1,jpj |
---|
| 268 | DO ji = 1,jpi |
---|
| 269 | d(1,1) = env_up(ji,jj) |
---|
| 270 | d(1,2) = env_dw(ji,jj) |
---|
| 271 | ! |
---|
| 272 | ! Computing derivative analytically |
---|
| 273 | ! |
---|
| 274 | ! 1. Derivative for upper envelope |
---|
| 275 | ibcbeg = 1 |
---|
| 276 | IF (nn_strt(je-1) == 2) THEN |
---|
| 277 | alpha = rn_e_al(je-1) |
---|
| 278 | num_s_up = nn_slev(je-1) |
---|
| 279 | coefa = rn_e_th(je-1) / & |
---|
| 280 | (env1(ji,jj)-env0(ji,jj)) |
---|
| 281 | coefb = rn_e_bb(je-1) + & |
---|
| 282 | ((env1(ji,jj)-env0(ji,jj))*rn_e_ba(je-1)) |
---|
| 283 | coefb = 1.0_wp-(coefb/(env1(ji,jj)-env0(ji,jj))) |
---|
| 284 | ELSE |
---|
| 285 | alpha = 0._wp |
---|
| 286 | num_s_up = 1 |
---|
| 287 | coefa = rn_e_th(je-1) |
---|
| 288 | coefb = rn_e_bb(je-1) |
---|
| 289 | END IF |
---|
| 290 | D1s_up = D1s_coord(-1._wp, coefa, coefb, & |
---|
| 291 | alpha, num_s_up, nn_strt(je-1)) |
---|
| 292 | IF (nn_strt(je-1) == 2) THEN |
---|
| 293 | d(2,1) = D1z_mes(env0(ji,jj), env1(ji,jj), D1s_up, & |
---|
| 294 | 0._wp, 1) |
---|
| 295 | ELSE |
---|
| 296 | d(2,1) = D1z_mes(env0(ji,jj), env1(ji,jj), D1s_up, & |
---|
| 297 | rn_e_hc(je-1), 1) |
---|
| 298 | END IF |
---|
| 299 | ! |
---|
| 300 | ! 2. Derivative for lower envelope |
---|
[15131] | 301 | IF ( je < tot_env ) THEN |
---|
[15126] | 302 | ibcend = 1 ! Bound. cond for the 1st derivative |
---|
| 303 | IF (nn_strt(je+1) == 2) THEN |
---|
| 304 | alpha = rn_e_al(je+1) |
---|
| 305 | num_s_dw = nn_slev(je+1) |
---|
| 306 | coefa = rn_e_th(je+1) / & |
---|
| 307 | (env3(ji,jj)-env2(ji,jj)) |
---|
| 308 | coefb = rn_e_bb(je+1) + & |
---|
| 309 | ((env3(ji,jj)-env2(ji,jj))*rn_e_ba(je+1)) |
---|
| 310 | coefb = 1.0_wp-(coefb/(env3(ji,jj)-env2(ji,jj))) |
---|
| 311 | ELSE |
---|
| 312 | alpha = 0._wp |
---|
| 313 | num_s_dw = 1 |
---|
| 314 | coefa = rn_e_th(je+1) |
---|
| 315 | coefb = rn_e_bb(je+1) |
---|
| 316 | END IF |
---|
| 317 | D1s_dw = D1s_coord(0._wp, coefa, coefb, & |
---|
| 318 | alpha, num_s_dw, nn_strt(je+1)) |
---|
| 319 | |
---|
| 320 | IF (nn_strt(je+1) == 2) THEN |
---|
| 321 | d(2,2) = D1z_mes(env2(ji,jj), env3(ji,jj), D1s_dw, & |
---|
| 322 | 0._wp, 1) |
---|
| 323 | ELSE |
---|
| 324 | d(2,2) = D1z_mes(env2(ji,jj), env3(ji,jj), D1s_dw, & |
---|
| 325 | rn_e_hc(je+1), 1) |
---|
| 326 | END IF |
---|
| 327 | ELSE |
---|
| 328 | ibcend = 2 ! Bound. cond for the 2nd derivative |
---|
| 329 | d(2,2) = 0._wp ! 2nd der = 0 |
---|
| 330 | ENDIF |
---|
| 331 | |
---|
| 332 | ! Computing spline's coefficients |
---|
| 333 | !IF ( lwp ) THEN |
---|
| 334 | ! WRITE(numout,*) "BOUNDARY COND. to CONSTRAIN SPLINES:" |
---|
| 335 | ! WRITE(numout,*) "ibcbeg: ", ibcbeg |
---|
| 336 | ! WRITE(numout,*) "ibcend: ", ibcend |
---|
| 337 | !ENDIF |
---|
[15134] | 338 | |
---|
[15126] | 339 | c = cub_spl(tau, d, n, ibcbeg, ibcend ) |
---|
| 340 | |
---|
[15134] | 341 | env_hc = rn_e_hc(je) |
---|
| 342 | coefa = rn_e_th(je) |
---|
| 343 | coefb = rn_e_bb(je) |
---|
| 344 | alpha = rn_e_al(je) |
---|
| 345 | |
---|
[15126] | 346 | DO jk = 1, num_s |
---|
| 347 | ind = s_1st+jk-1 |
---|
[15134] | 348 | ! |
---|
| 349 | ! Computing stretched distribution of interpolation points |
---|
| 350 | ! W-GRID |
---|
| 351 | zcoefw = -s_coord( sigma(1,jk,num_s), coefa, & |
---|
| 352 | coefb, alpha, num_s, nn_strt(je) ) |
---|
| 353 | ! T-GRID |
---|
| 354 | zcoeft = -s_coord( sigma(0,jk,num_s), coefa, & |
---|
| 355 | coefb, alpha, num_s, nn_strt(je) ) |
---|
| 356 | ! |
---|
| 357 | z_gsigw3(ji,jj,ind) = zcoefw |
---|
| 358 | z_gsigt3(ji,jj,ind) = zcoeft |
---|
| 359 | ! |
---|
| 360 | gdept_0(ji,jj,ind) = z_cub_spl(zcoeft, tau, c) |
---|
| 361 | gdepw_0(ji,jj,ind) = z_cub_spl(zcoefw, tau, c) |
---|
[15126] | 362 | END DO |
---|
| 363 | |
---|
| 364 | END DO |
---|
| 365 | END DO |
---|
| 366 | |
---|
| 367 | END IF |
---|
| 368 | |
---|
| 369 | END DO |
---|
| 370 | |
---|
| 371 | ! ============================================= |
---|
| 372 | ! 4. COMPUTE e3t, e3w for ALL general levels |
---|
| 373 | ! as finite differences |
---|
| 374 | ! ============================================= |
---|
| 375 | ! |
---|
| 376 | DO jk=1,jpkm1 |
---|
| 377 | e3t_0(:,:,jk) = gdepw_0(:,:,jk+1) - gdepw_0(:,:,jk) |
---|
| 378 | e3w_0(:,:,jk+1) = gdept_0(:,:,jk+1) - gdept_0(:,:,jk) |
---|
| 379 | ENDDO |
---|
| 380 | ! Surface |
---|
| 381 | jk = 1 |
---|
| 382 | e3w_0(:,:,jk) = 2.0_wp * (gdept_0(:,:,1) - gdepw_0(:,:,1)) |
---|
| 383 | ! |
---|
| 384 | ! Bottom |
---|
| 385 | jk = jpk |
---|
| 386 | e3t_0(:,:,jk) = 2.0_wp * (gdept_0(:,:,jk) - gdepw_0(:,:,jk)) |
---|
| 387 | |
---|
| 388 | !============================== |
---|
| 389 | ! Computing mbathy |
---|
| 390 | !============================== |
---|
| 391 | |
---|
| 392 | IF( lwp ) WRITE(numout,*) '' |
---|
| 393 | IF( lwp ) WRITE(numout,*) 'MES mbathy:' |
---|
| 394 | IF( lwp ) WRITE(numout,*) '' |
---|
| 395 | |
---|
| 396 | DO jj = 1, jpj |
---|
| 397 | DO ji = 1, jpi |
---|
| 398 | DO jk = 1, jpkm1 |
---|
| 399 | !IF( lwp ) WRITE(numout,*) 'scobot: ', scobot(ji,jj), 'gdept_0: ', gdept_0(ji,jj,jk) |
---|
| 400 | IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) |
---|
| 401 | IF( scobot(ji,jj) == 0.e0 ) mbathy(ji,jj) = 0 |
---|
| 402 | IF( scobot(ji,jj) < 0.e0 ) mbathy(ji,jj) = int( scobot(ji,jj)) !???? do we need it? |
---|
| 403 | END DO |
---|
| 404 | !IF( lwp ) WRITE(numout,*) 'mbathy(',ji,',',jj,') = ', mbathy(ji,jj) |
---|
| 405 | END DO |
---|
| 406 | END DO |
---|
| 407 | |
---|
| 408 | !============================================== |
---|
| 409 | ! Computing e3u_0, e3v_0, e3f_0, e3uw_0, e3vw_0 |
---|
| 410 | !============================================== |
---|
| 411 | |
---|
| 412 | DO jj = 1, jpjm1 |
---|
| 413 | DO ji = 1, jpim1 |
---|
| 414 | DO jk = 1, jpk |
---|
| 415 | |
---|
| 416 | e3u_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3t_0(ji,jj,jk) + & |
---|
| 417 | MIN(1,mbathy(ji+1,jj))*e3t_0(ji+1,jj,jk) ) / & |
---|
| 418 | MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji+1,jj)) ) |
---|
| 419 | |
---|
| 420 | e3v_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3t_0(ji,jj,jk) + & |
---|
| 421 | MIN(1,mbathy(ji,jj+1))*e3t_0(ji,jj+1,jk) ) / & |
---|
| 422 | MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1)) ) |
---|
| 423 | |
---|
| 424 | e3uw_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3w_0(ji,jj,jk) + & |
---|
| 425 | MIN(1,mbathy(ji+1,jj))*e3w_0(ji+1,jj,jk) ) / & |
---|
| 426 | MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji+1,jj)) ) |
---|
| 427 | |
---|
| 428 | e3vw_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3w_0(ji,jj,jk) + & |
---|
| 429 | MIN(1,mbathy(ji,jj+1))*e3w_0(ji,jj+1,jk) ) / & |
---|
| 430 | MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1)) ) |
---|
| 431 | |
---|
| 432 | e3f_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3t_0(ji,jj,jk) + & |
---|
| 433 | MIN(1,mbathy(ji+1,jj))*e3t_0(ji+1,jj,jk) + & |
---|
| 434 | MIN(1,mbathy(ji+1,jj+1))* e3t_0(ji+1,jj+1,jk) + & |
---|
| 435 | MIN(1,mbathy(ji,jj+1))*e3t_0(ji,jj+1,jk) ) / & |
---|
| 436 | MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1)) & |
---|
| 437 | + MIN(1,mbathy(ji+1,jj))+MIN(1,mbathy(ji+1,jj+1)) ) |
---|
| 438 | ENDDO |
---|
| 439 | ENDDO |
---|
| 440 | ENDDO |
---|
| 441 | |
---|
| 442 | ! Adjusting for geopotential levels, if any |
---|
| 443 | |
---|
| 444 | DO je = 1, tot_env ! LOOP over all the requested envelopes |
---|
| 445 | |
---|
| 446 | IF (je == 1) THEN |
---|
| 447 | s_1st = 1 |
---|
| 448 | num_s = nn_slev(je) |
---|
| 449 | max_env_up = 0._wp ! surface |
---|
| 450 | min_env_up = 0._wp ! surface |
---|
| 451 | ELSE |
---|
| 452 | s_1st = s_1st + num_s - 1 |
---|
| 453 | num_s = nn_slev(je) + 1 |
---|
[15131] | 454 | max_env_up = max_dep_env(je-1) |
---|
| 455 | min_env_up = min_dep_env(je-1) |
---|
[15126] | 456 | ENDIF |
---|
| 457 | |
---|
[15131] | 458 | max_env_dw = max_dep_env(je) |
---|
| 459 | min_env_dw = min_dep_env(je) |
---|
[15126] | 460 | |
---|
| 461 | IF ( max_env_up == min_env_up .AND. max_env_dw == min_env_dw ) THEN |
---|
| 462 | |
---|
| 463 | DO jj = 1,jpjm1 |
---|
| 464 | DO ji = 1,jpim1 |
---|
| 465 | DO jk = 1, num_s |
---|
| 466 | |
---|
| 467 | ind = s_1st+jk-1 |
---|
| 468 | e3u_0 (ji,jj,ind) = MIN( e3t_0(ji,jj,ind), e3t_0(ji+1,jj,ind)) |
---|
| 469 | e3uw_0(ji,jj,ind) = MIN( e3w_0(ji,jj,ind), e3w_0(ji+1,jj,ind)) |
---|
| 470 | e3v_0 (ji,jj,ind) = MIN( e3t_0(ji,jj,ind), e3t_0(ji,jj+1,ind)) |
---|
| 471 | e3vw_0(ji,jj,ind) = MIN( e3w_0(ji,jj,ind), e3w_0(ji,jj+1,ind)) |
---|
| 472 | e3f_0 (ji,jj,ind) = MIN( e3t_0(ji,jj,ind), e3t_0(ji+1,jj,ind)) |
---|
| 473 | |
---|
| 474 | ENDDO |
---|
| 475 | ENDDO |
---|
| 476 | ENDDO |
---|
| 477 | |
---|
| 478 | ENDIF |
---|
| 479 | |
---|
| 480 | ENDDO |
---|
| 481 | |
---|
| 482 | !============================== |
---|
| 483 | ! Computing gdept3w |
---|
| 484 | !============================== |
---|
| 485 | |
---|
| 486 | gde3w_0(:,:,1) = 0.5 * e3w_0(:,:,1) |
---|
| 487 | DO jk = 2, jpk |
---|
| 488 | gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk) |
---|
| 489 | END DO |
---|
| 490 | |
---|
| 491 | !================================================================================= |
---|
| 492 | ! From here equal to sco code - domzgr.F90 line 2079 |
---|
| 493 | ! Lateral B.C. |
---|
| 494 | |
---|
| 495 | CALL lbc_lnk( e3t_0 , 'T', 1._wp ) |
---|
| 496 | CALL lbc_lnk( e3u_0 , 'U', 1._wp ) |
---|
| 497 | CALL lbc_lnk( e3v_0 , 'V', 1._wp ) |
---|
| 498 | CALL lbc_lnk( e3f_0 , 'F', 1._wp ) |
---|
| 499 | CALL lbc_lnk( e3w_0 , 'W', 1._wp ) |
---|
| 500 | CALL lbc_lnk( e3uw_0, 'U', 1._wp ) |
---|
| 501 | CALL lbc_lnk( e3vw_0, 'V', 1._wp ) |
---|
| 502 | |
---|
| 503 | where (e3t_0 (:,:,:).eq.0.0) e3t_0(:,:,:) = 1.0 |
---|
| 504 | where (e3u_0 (:,:,:).eq.0.0) e3u_0(:,:,:) = 1.0 |
---|
| 505 | where (e3v_0 (:,:,:).eq.0.0) e3v_0(:,:,:) = 1.0 |
---|
| 506 | where (e3f_0 (:,:,:).eq.0.0) e3f_0(:,:,:) = 1.0 |
---|
| 507 | where (e3w_0 (:,:,:).eq.0.0) e3w_0(:,:,:) = 1.0 |
---|
| 508 | where (e3uw_0 (:,:,:).eq.0.0) e3uw_0(:,:,:) = 1.0 |
---|
| 509 | where (e3vw_0 (:,:,:).eq.0.0) e3vw_0(:,:,:) = 1.0 |
---|
| 510 | |
---|
| 511 | IF( lwp ) WRITE(numout,*) 'Refine mbathy' |
---|
| 512 | DO jj = 1, jpj |
---|
| 513 | DO ji = 1, jpi |
---|
| 514 | DO jk = 1, jpkm1 |
---|
| 515 | IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) |
---|
| 516 | END DO |
---|
| 517 | IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0 |
---|
| 518 | END DO |
---|
| 519 | END DO |
---|
| 520 | |
---|
| 521 | IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), & |
---|
| 522 | & ' MAX ', MAXVAL( mbathy(:,:) ) |
---|
| 523 | |
---|
| 524 | IF( nprint == 1 .AND. lwp ) THEN ! min max values over the local domain |
---|
| 525 | WRITE(numout,*) 'zgr_sco : mbathy' |
---|
| 526 | WRITE(numout,*) '~~~~~~~' |
---|
| 527 | WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) |
---|
| 528 | WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), & |
---|
| 529 | & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ' , MINVAL( gde3w_0(:,:,:) ) |
---|
| 530 | WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 (:,:,:) ), & |
---|
| 531 | & ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 (:,:,:) ), & |
---|
| 532 | & ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 (:,:,:) ), & |
---|
| 533 | & ' w ', MINVAL( e3w_0 (:,:,:) ) |
---|
| 534 | |
---|
| 535 | WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), & |
---|
| 536 | & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ' , MAXVAL( gde3w_0(:,:,:) ) |
---|
| 537 | WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 (:,:,:) ), & |
---|
| 538 | & ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 (:,:,:) ), & |
---|
| 539 | & ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 (:,:,:) ), & |
---|
| 540 | & ' w ', MAXVAL( e3w_0 (:,:,:) ) |
---|
| 541 | ENDIF |
---|
| 542 | |
---|
| 543 | !================================================================================ |
---|
| 544 | ! Check coordinates makes sense |
---|
| 545 | !================================================================================ |
---|
| 546 | |
---|
[15134] | 547 | ! Extracting MEs depth profile in the shallowest point of the deepest |
---|
| 548 | ! envelope for a first check of monotonicity of transformation |
---|
| 549 | ! (also useful for debugging) |
---|
| 550 | pmsk(:,:) = 0.0 |
---|
| 551 | WHERE ( bathy > 0.0 ) pmsk = 1.0 |
---|
| 552 | CALL mpp_minloc( envlt(:,:,tot_env), pmsk(:,:), mindep, eii, eij ) |
---|
[15172] | 553 | IF ((mi0(eii)>=1 .AND. mi0(eii)<=jpi) .AND. (mj0(eij)>=1 .AND. mj0(eij)<=jpj)) THEN |
---|
[15134] | 554 | irnk = mpprank |
---|
| 555 | DO je = 1, tot_env |
---|
| 556 | rn_ebot_max(je) = envlt(mi0(eii),mj0(eij),je) |
---|
| 557 | END DO |
---|
| 558 | z_gsigt1(:) = z_gsigt3(mi0(eii),mj0(eij),:) |
---|
| 559 | z_gsigw1(:) = z_gsigw3(mi0(eii),mj0(eij),:) |
---|
| 560 | gdept1(:) = gdept_0(mi0(eii),mj0(eij),:) |
---|
| 561 | gdepw1(:) = gdepw_0(mi0(eii),mj0(eij),:) |
---|
| 562 | e3t1(:) = e3t_0(mi0(eii),mj0(eij),:) |
---|
| 563 | e3w1(:) = e3w_0(mi0(eii),mj0(eij),:) |
---|
| 564 | ELSE |
---|
| 565 | irnk = -1 |
---|
| 566 | END IF |
---|
| 567 | IF( lk_mpp ) CALL mppsync |
---|
| 568 | IF( lk_mpp ) CALL mpp_max(irnk) |
---|
| 569 | IF( lk_mpp ) CALL mppbcast_a_real(rn_ebot_max, max_nn_env, irnk ) |
---|
| 570 | IF( lk_mpp ) CALL mppbcast_a_real(z_gsigt1 , jpk , irnk ) |
---|
| 571 | IF( lk_mpp ) CALL mppbcast_a_real(z_gsigw1 , jpk , irnk ) |
---|
| 572 | IF( lk_mpp ) CALL mppbcast_a_real(gdept1 , jpk , irnk ) |
---|
| 573 | IF( lk_mpp ) CALL mppbcast_a_real(gdepw1 , jpk , irnk ) |
---|
| 574 | IF( lk_mpp ) CALL mppbcast_a_real(e3t1 , jpk , irnk ) |
---|
| 575 | IF( lk_mpp ) CALL mppbcast_a_real(e3w1 , jpk , irnk ) |
---|
| 576 | |
---|
| 577 | IF( lwp ) THEN |
---|
| 578 | WRITE(numout,*) "" |
---|
| 579 | WRITE(numout,*) "mes_build:" |
---|
| 580 | WRITE(numout,*) "~~~~~~~~~" |
---|
| 581 | WRITE(numout,*) "" |
---|
| 582 | WRITE(numout,*) " FIRST CHECK: Checking MEs-levels profile in the shallowest point of the last envelope:" |
---|
| 583 | WRITE(numout,*) " it is the most likely point where monotonicty of splines may be violeted. " |
---|
| 584 | WRITE(numout,*) "" |
---|
| 585 | DO je = 1, tot_env |
---|
| 586 | WRITE(numout,*) ' * depth of envelope ', je, ' at point (',eii,',',eij,') is ', rn_ebot_max(je) |
---|
| 587 | END DO |
---|
| 588 | WRITE(numout,*) "" |
---|
| 589 | WRITE(numout,*) " MEs-coordinates" |
---|
| 590 | WRITE(numout,*) "" |
---|
| 591 | WRITE(numout,*) " k z_gsigw1 z_gsigt1" |
---|
| 592 | WRITE(numout,*) "" |
---|
| 593 | DO jk = 1, jpk |
---|
| 594 | WRITE(numout,*) ' ', jk, z_gsigw1(jk), z_gsigt1(jk) |
---|
| 595 | ENDDO |
---|
| 596 | WRITE(numout,*) "" |
---|
| 597 | WRITE(numout,*) "-----------------------------------------------------------------" |
---|
| 598 | WRITE(numout,*) "" |
---|
| 599 | WRITE(numout,*) " MEs-levels depths and scale factors" |
---|
| 600 | WRITE(numout,*) "" |
---|
| 601 | WRITE(numout,*) " k gdepw1 e3w1 gdept1 e3t1" |
---|
| 602 | WRITE(numout,*) "" |
---|
| 603 | DO jk = 1, jpk |
---|
| 604 | WRITE(numout,*) ' ', jk, gdepw1(jk), e3w1(jk), gdept1(jk), e3t1(jk) |
---|
| 605 | ENDDO |
---|
| 606 | WRITE(numout,*) "-----------------------------------------------------------------" |
---|
| 607 | ENDIF |
---|
| 608 | |
---|
| 609 | ! Checking monotonicity for cubic splines |
---|
| 610 | DO jk = 1, jpk-1 |
---|
| 611 | IF ( gdept1(jk+1) < gdept1(jk) ) THEN |
---|
| 612 | WRITE(ctlmes,*) 'NOT MONOTONIC gdept_0: change envelopes' |
---|
| 613 | CALL ctl_stop( ctlmes ) |
---|
| 614 | END IF |
---|
| 615 | IF ( gdepw1(jk+1) < gdepw1(jk) ) THEN |
---|
| 616 | WRITE(ctlmes,*) 'NOT MONOTONIC gdepw_0: change envelopes' |
---|
| 617 | CALL ctl_stop( ctlmes ) |
---|
| 618 | END IF |
---|
| 619 | IF ( e3t1(jk) < 0.0 ) THEN |
---|
| 620 | WRITE(ctlmes,*) 'NEGATIVE e3t_0: change envelopes' |
---|
| 621 | CALL ctl_stop( ctlmes ) |
---|
| 622 | END IF |
---|
| 623 | IF ( e3w1(jk) < 0.0 ) THEN |
---|
| 624 | WRITE(ctlmes,*) 'NEGATIVE e3w_0: change envelopes' |
---|
| 625 | CALL ctl_stop( ctlmes ) |
---|
| 626 | END IF |
---|
| 627 | END DO |
---|
| 628 | |
---|
| 629 | ! CHECKING THE WHOLE DOMAIN |
---|
[15126] | 630 | DO ji = 1, jpi |
---|
| 631 | DO jj = 1, jpj |
---|
[15134] | 632 | DO jk = 1, jpk !mbathy(ji,jj) |
---|
[15126] | 633 | ! check coordinate is monotonically increasing |
---|
| 634 | IF (e3w_0(ji,jj,jk) <= 0._wp .OR. e3t_0(ji,jj,jk) <= 0._wp ) THEN |
---|
[15134] | 635 | WRITE(ctmp1,*) 'ERROR mes_build: e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk |
---|
| 636 | WRITE(numout,*) 'ERROR mes_build: e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk |
---|
[15126] | 637 | WRITE(numout,*) 'e3w',e3w_0(ji,jj,:) |
---|
| 638 | WRITE(numout,*) 'e3t',e3t_0(ji,jj,:) |
---|
| 639 | CALL ctl_stop( ctmp1 ) |
---|
| 640 | ENDIF |
---|
| 641 | ! and check it has never gone negative |
---|
| 642 | IF ( gdepw_0(ji,jj,jk) < 0._wp .OR. gdept_0(ji,jj,jk) < 0._wp ) THEN |
---|
[15134] | 643 | WRITE(ctmp1,*) 'ERROR mes_build: gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk |
---|
| 644 | WRITE(numout,*) 'ERROR mes_build: gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk |
---|
[15126] | 645 | WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:) |
---|
| 646 | WRITE(numout,*) 'gdept',gdept_0(ji,jj,:) |
---|
| 647 | CALL ctl_stop( ctmp1 ) |
---|
| 648 | ENDIF |
---|
| 649 | END DO |
---|
| 650 | END DO |
---|
| 651 | END DO |
---|
| 652 | ! |
---|
| 653 | CALL wrk_dealloc( jpi, jpj, max_nn_env , envlt ) |
---|
| 654 | ! |
---|
[15129] | 655 | IF( nn_timing == 1 ) CALL timing_stop('mes_build') |
---|
[15126] | 656 | |
---|
[15129] | 657 | END SUBROUTINE mes_build |
---|
[15126] | 658 | |
---|
| 659 | ! ===================================================================================================== |
---|
| 660 | |
---|
[15129] | 661 | SUBROUTINE mes_init |
---|
[15134] | 662 | !!----------------------------------------------------------------------------- |
---|
| 663 | !! *** ROUTINE mes_init *** |
---|
| 664 | !! |
---|
| 665 | !! ** Purpose : Initilaise a Multi Enveloped s-coordinate (MEs) system |
---|
| 666 | !! |
---|
| 667 | !!----------------------------------------------------------------------------- |
---|
| 668 | |
---|
[15129] | 669 | CHARACTER(lc) :: env_name ! name of the externally defined envelope |
---|
[15134] | 670 | INTEGER :: ji, jj, je, inum |
---|
| 671 | INTEGER :: cor_env, ios |
---|
| 672 | REAL(wp) :: rn_bot_min ! min depth of ocean bottom (>0) (m) |
---|
| 673 | REAL(wp) :: rn_bot_max ! max depth of ocean bottom (= ocean depth) (>0) (m) |
---|
[15129] | 674 | |
---|
[15278] | 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 |
---|
[15129] | 680 | !!---------------------------------------------------------------------- |
---|
| 681 | ! |
---|
| 682 | IF( nn_timing == 1 ) CALL timing_start('mes_init') |
---|
| 683 | ! |
---|
| 684 | CALL wrk_alloc( jpi, jpj, max_nn_env , envlt ) |
---|
| 685 | ! |
---|
[15131] | 686 | ! Initialising some arrays |
---|
| 687 | max_dep_env(:) = 0.0 |
---|
| 688 | min_dep_env(:) = 0.0 |
---|
| 689 | ! |
---|
| 690 | IF(lwp) THEN ! control print |
---|
| 691 | WRITE(numout,*) '' |
---|
| 692 | WRITE(numout,*) 'mes_init: Setting a Multi-Envelope s-ccordinate system (Bruciaferri et al. 2018)' |
---|
| 693 | WRITE(numout,*) '~~~~~~~~' |
---|
| 694 | END IF |
---|
| 695 | ! |
---|
[15129] | 696 | ! Namelist namzgr_mes in reference namelist : envelopes and |
---|
| 697 | ! sigma-stretching parameters |
---|
| 698 | REWIND( numnam_ref ) |
---|
| 699 | READ ( numnam_ref, namzgr_mes, IOSTAT = ios, ERR = 901) |
---|
| 700 | 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_mes in reference namelist', lwp ) |
---|
| 701 | ! Namelist namzgr_mes in configuration namelist : envelopes and |
---|
| 702 | ! sigma-stretching parameters |
---|
| 703 | REWIND( numnam_cfg ) |
---|
| 704 | READ ( numnam_cfg, namzgr_mes, IOSTAT = ios, ERR = 902 ) |
---|
| 705 | 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_mes in configuration namelist', lwp ) |
---|
| 706 | IF(lwm) WRITE ( numond, namzgr_mes ) |
---|
| 707 | ! |
---|
| 708 | ! 1) Check if namelist is defined correctly. |
---|
| 709 | ! Not strictly needed but we force the user |
---|
| 710 | ! to define the namelist correctly. |
---|
| 711 | tot_env = 0 ! total number of requested envelopes |
---|
| 712 | cor_env = 0 ! total number of correctly defined envelopes |
---|
| 713 | DO je = 1, max_nn_env |
---|
| 714 | IF ( ln_envl(je) ) tot_env = tot_env + 1 |
---|
| 715 | IF ( ln_envl(je) .AND. cor_env == (je-1)) cor_env = cor_env + 1 |
---|
| 716 | ENDDO |
---|
| 717 | WRITE(ctlmes,*) 'num. of REQUESTED env. and num. of CORRECTLY defined env. are DIFFERENT' |
---|
| 718 | IF ( tot_env /= cor_env ) CALL ctl_stop( ctlmes ) |
---|
[15131] | 719 | IF(lwp) THEN |
---|
| 720 | WRITE(numout,*) ' Num. or requested envelopes: ', tot_env |
---|
| 721 | WRITE(numout,*) '' |
---|
| 722 | END IF |
---|
[15129] | 723 | ! |
---|
| 724 | ! 2) Checking consistency of user defined parameters |
---|
| 725 | WRITE(ctlmes,*) 'TOT number of levels (jpk) IS DIFFERENT from sum over nn_slev(:)' |
---|
| 726 | IF ( SUM(nn_slev(:)) /= jpk ) CALL ctl_stop( ctlmes ) |
---|
| 727 | ! |
---|
[15278] | 728 | IF( .NOT.(ln_loc_zgr) .AND. ln_pst_l2g) ln_pst_l2g = .FALSE. |
---|
| 729 | ! |
---|
[15129] | 730 | ! 3) Reading Bathymetry and envelopes |
---|
| 731 | IF( ntopo == 1 ) THEN |
---|
[15131] | 732 | IF(lwp) THEN |
---|
| 733 | WRITE(numout,*) ' Reading bathymetry and envelopes ...' |
---|
| 734 | WRITE(numout,*) '' |
---|
| 735 | END IF |
---|
| 736 | |
---|
[15129] | 737 | CALL iom_open ( 'bathy_meter.nc', inum ) |
---|
| 738 | CALL iom_get ( inum, jpdom_data, 'Bathymetry' , bathy, lrowattr=ln_use_jattr ) |
---|
| 739 | DO je = 1, tot_env |
---|
| 740 | WRITE (env_name, '(A6,I0)') 'hbatt_', je |
---|
| 741 | CALL iom_get ( inum, jpdom_data, TRIM(env_name) , envlt(:,:,je), lrowattr=ln_use_jattr ) |
---|
| 742 | ENDDO |
---|
| 743 | CALL iom_close( inum ) |
---|
| 744 | ELSE |
---|
| 745 | WRITE(ctlmes,*) 'parameter , ntopo = ', ntopo |
---|
| 746 | CALL ctl_stop( ctlmes ) |
---|
| 747 | ENDIF |
---|
| 748 | IF( nn_closea == 0 ) CALL clo_bat( bathy, mbathy ) !== NO closed seas or lakes ==! |
---|
| 749 | ! |
---|
| 750 | ! 4) Checking consistency of envelopes |
---|
| 751 | DO je = 1, tot_env-1 |
---|
| 752 | WRITE(ctlmes,*) 'Envelope ', je+1, ' is shallower that Envelope ', je |
---|
| 753 | IF (MAXVAL(envlt(:,:,je+1)) < MAXVAL(envlt(:,:,je))) CALL ctl_stop( ctlmes ) |
---|
| 754 | ENDDO |
---|
[15134] | 755 | ! 5) Computing max and min depths of envelopes |
---|
[15129] | 756 | DO je = 1, tot_env |
---|
[15131] | 757 | max_dep_env(je) = MAXVAL(envlt(:,:,je)) |
---|
| 758 | min_dep_env(je) = MINVAL(envlt(:,:,je)) |
---|
| 759 | IF( lk_mpp ) CALL mpp_max( max_dep_env(je) ) |
---|
| 760 | IF( lk_mpp ) CALL mpp_min( min_dep_env(je) ) |
---|
[15129] | 761 | END DO |
---|
[15131] | 762 | IF( lk_mpp ) CALL mppsync |
---|
[15129] | 763 | ! |
---|
| 764 | ! 6) Set maximum and minimum ocean depth |
---|
| 765 | bathy(:,:) = MIN( rn_bot_max, bathy(:,:) ) |
---|
| 766 | DO jj = 1, jpj |
---|
| 767 | DO ji = 1, jpi |
---|
| 768 | IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_bot_min, bathy(ji,jj) ) |
---|
| 769 | END DO |
---|
| 770 | END DO |
---|
| 771 | ! |
---|
| 772 | IF(lwp) THEN ! control print |
---|
| 773 | WRITE(numout,*) |
---|
[15131] | 774 | WRITE(numout,*) ' GEOMETRICAL FEATURES OF THE REQUESTED MEs GRID:' |
---|
| 775 | WRITE(numout,*) ' -----------------------------------------------' |
---|
| 776 | WRITE(numout,*) ' Minimum depth of the ocean rn_bot_min = ', rn_bot_min |
---|
| 777 | WRITE(numout,*) ' Maximum depth of the ocean rn_bot_max = ', rn_bot_max |
---|
[15129] | 778 | WRITE(numout,*) '' |
---|
[15131] | 779 | WRITE(numout,*) ' ENVELOPES:' |
---|
| 780 | WRITE(numout,*) ' ========= ' |
---|
[15129] | 781 | WRITE(numout,*) '' |
---|
[15131] | 782 | DO je = 1, tot_env |
---|
| 783 | WRITE(numout,*) ' * envelope ', je, ': min depth = ', min_dep_env(je) |
---|
| 784 | WRITE(numout,*) ' max depth = ', max_dep_env(je) |
---|
| 785 | WRITE(numout,*) '' |
---|
| 786 | END DO |
---|
| 787 | WRITE(numout,*) '' |
---|
| 788 | WRITE(numout,*) ' SUBDOMAINS:' |
---|
| 789 | WRITE(numout,*) ' ========== ' |
---|
| 790 | WRITE(numout,*) '' |
---|
| 791 | DO je = 1, tot_env |
---|
| 792 | WRITE(numout,*) ' * subdomain ', je,':' |
---|
[15129] | 793 | IF ( je == 1) THEN |
---|
[15131] | 794 | WRITE(numout,*) ' Shallower envelope: free surface' |
---|
[15129] | 795 | ELSE |
---|
[15131] | 796 | WRITE(numout,*) ' Shallower envelope: envelope ',je-1 |
---|
[15129] | 797 | END IF |
---|
[15131] | 798 | WRITE(numout,*) ' Deeper envelope: envelope ',je |
---|
| 799 | WRITE(numout,*) ' Number of MEs-levs: nn_slev(',je,') = ',nn_slev(je) |
---|
[15129] | 800 | IF ( isodd(je) ) THEN |
---|
[15131] | 801 | WRITE(numout,*) ' Stretched s-coordinates: ' |
---|
[15129] | 802 | ELSE |
---|
[15131] | 803 | WRITE(numout,*) ' Stretched CUBIC SPLINES: ' |
---|
[15129] | 804 | END IF |
---|
[15131] | 805 | IF (nn_strt(je) == 0) WRITE(numout,*) ' M96 stretching function' |
---|
| 806 | IF (nn_strt(je) == 1) WRITE(numout,*) ' SH94 stretching function' |
---|
| 807 | IF (nn_strt(je) == 2) WRITE(numout,*) ' SF12 stretching function' |
---|
| 808 | WRITE(numout,*) ' critical depth rn_e_hc(',je,') = ',rn_e_hc(je) |
---|
| 809 | WRITE(numout,*) ' surface stretc. coef. rn_e_th(',je,') = ',rn_e_th(je) |
---|
[15129] | 810 | IF (nn_strt(je) == 2) THEN |
---|
[15131] | 811 | WRITE(numout,*) ' bottom stretc. coef. rn_e_ba(',je,') = ',rn_e_ba(je) |
---|
[15129] | 812 | END IF |
---|
[15131] | 813 | WRITE(numout,*) ' bottom stretc. coef. rn_e_bb(',je,') = ',rn_e_bb(je) |
---|
[15129] | 814 | IF (nn_strt(je) == 2) THEN |
---|
[15131] | 815 | WRITE(numout,*) ' bottom stretc. coef. rn_e_al(',je,') = ',rn_e_al(je) |
---|
[15129] | 816 | END IF |
---|
[15131] | 817 | WRITE(numout,*) ' ------------------------------------------------------------------' |
---|
[15129] | 818 | ENDDO |
---|
[15278] | 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 |
---|
[15129] | 840 | ENDIF |
---|
| 841 | |
---|
| 842 | END SUBROUTINE mes_init |
---|
| 843 | |
---|
| 844 | ! ===================================================================================================== |
---|
| 845 | |
---|
[15126] | 846 | FUNCTION isodd( a ) RESULT(odd) |
---|
| 847 | ! Notice that this method uses a btest intrinsic function. If you look at any number in binary |
---|
| 848 | ! notation, you'll note that the Least Significant Bit (LSB) will always be set for odd numbers |
---|
| 849 | ! and cleared for even numbers. Hence, all that is necessary to determine if a number is even or |
---|
| 850 | ! odd is to check the LSB. Since bitwise operations like AND, OR, XOR etc. are usually much faster |
---|
| 851 | ! than the arithmetic operations, this method will run a lot faster than the one with modulo |
---|
| 852 | ! operation. |
---|
| 853 | ! http://forums.devshed.com/software-design-43/quick-algorithm-determine-odd-29843.html |
---|
| 854 | |
---|
| 855 | INTEGER, INTENT (in) :: a |
---|
| 856 | LOGICAL :: odd |
---|
| 857 | odd = btest(a, 0) |
---|
| 858 | END FUNCTION isodd |
---|
| 859 | |
---|
| 860 | ! ===================================================================================================== |
---|
| 861 | |
---|
| 862 | FUNCTION iseven( a ) RESULT(even) |
---|
| 863 | INTEGER, INTENT (in) :: a |
---|
| 864 | LOGICAL :: even |
---|
| 865 | even = .NOT. isodd(a) |
---|
| 866 | END FUNCTION iseven |
---|
| 867 | |
---|
| 868 | ! ===================================================================================================== |
---|
| 869 | |
---|
| 870 | FUNCTION sech( x ) RESULT( seh ) |
---|
| 871 | |
---|
| 872 | REAL(wp), INTENT(in ) :: x |
---|
| 873 | REAL(wp) :: seh |
---|
| 874 | |
---|
| 875 | seh = 1._wp / COSH(x) |
---|
| 876 | |
---|
| 877 | END FUNCTION sech |
---|
| 878 | |
---|
| 879 | ! ===================================================================================================== |
---|
| 880 | |
---|
| 881 | FUNCTION sigma( vgrid, pk, kmax ) RESULT( ps1 ) |
---|
| 882 | !!---------------------------------------------------------------------- |
---|
| 883 | !! *** ROUTINE sigma *** |
---|
| 884 | !! |
---|
| 885 | !! ** Purpose : provide the analytical function for sigma-coordinate |
---|
| 886 | !! (not stretched s-coordinate). |
---|
| 887 | !! |
---|
| 888 | !! ** Method : the function provide the non-dimensional position of |
---|
| 889 | !! T and W points (i.e. between 0 and 1). |
---|
| 890 | !! T-points at integer values (between 1 and jpk) |
---|
| 891 | !! W-points at integer values - 1/2 (between 0.5 and |
---|
| 892 | !! jpk-0.5) |
---|
| 893 | !! |
---|
| 894 | !!---------------------------------------------------------------------- |
---|
| 895 | INTEGER, INTENT (in) :: pk ! continuous k coordinate |
---|
| 896 | INTEGER, INTENT (in) :: kmax ! number of levels |
---|
| 897 | INTEGER, INTENT (in) :: vgrid ! type of vertical grid: 0 = T, 1 = W |
---|
| 898 | REAL(wp) :: kindx ! index of T or W points |
---|
| 899 | REAL(wp) :: ps1 ! value of sigma coordinate (-1 <= ps1 <=0) |
---|
| 900 | ! |
---|
| 901 | SELECT CASE (vgrid) |
---|
| 902 | |
---|
| 903 | CASE (0) ! T-points at integer values (between 1 and jpk) |
---|
| 904 | kindx = REAL(pk,wp) |
---|
| 905 | CASE (1) ! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) |
---|
| 906 | kindx = REAL(pk,wp) - 0.5_wp |
---|
| 907 | CASE DEFAULT |
---|
| 908 | WRITE(ctlmes,*) 'No valid vertical grid option selected' |
---|
| 909 | |
---|
| 910 | END SELECT |
---|
| 911 | |
---|
| 912 | |
---|
| 913 | ps1 = -(kindx - 0.5) / REAL(kmax-1) ! The surface is at the first W level |
---|
| 914 | ! while the bottom coincides with the |
---|
| 915 | ! last W level |
---|
| 916 | END FUNCTION sigma |
---|
| 917 | |
---|
| 918 | ! ===================================================================================================== |
---|
| 919 | |
---|
| 920 | FUNCTION s_coord( s, ca, cb, alpha, kmax, ftype ) RESULT ( pf1 ) |
---|
| 921 | !!---------------------------------------------------------------------- |
---|
| 922 | !! *** ROUTINE s_coord *** |
---|
| 923 | !! |
---|
| 924 | !! ** Purpose : provide the analytical stretching function |
---|
| 925 | !! for s-coordinate. |
---|
| 926 | !! |
---|
| 927 | !! ** Method : if ftype = 2: Siddorn and Furner 2012 |
---|
| 928 | !! analytical function is used |
---|
| 929 | !! if ftype = 1: Song and Haidvogel 1994 |
---|
| 930 | !! analytical function is used |
---|
| 931 | !! if ftype = 0: Madec et al. 1996 |
---|
| 932 | !! analytical function is used |
---|
| 933 | !! (pag 65 of NEMO Manual) |
---|
| 934 | !! Reference : Madec, Lott, Delecluse and |
---|
| 935 | !! Crepon, 1996. JPO, 26, 1393-1408 |
---|
| 936 | !! |
---|
| 937 | !! s MUST be NEGATIVE: -1 <= s <= 0 |
---|
| 938 | !!---------------------------------------------------------------------- |
---|
| 939 | REAL(wp), INTENT(in) :: s ! continuous not stretched |
---|
| 940 | ! "s" coordinate |
---|
| 941 | REAL(wp), INTENT(in) :: ca ! surface stretch. coeff |
---|
| 942 | REAL(wp), INTENT(in) :: cb ! bottom stretch. coeff |
---|
| 943 | REAL(wp), INTENT(in) :: alpha ! alpha stretch. coeff for SF12 |
---|
| 944 | INTEGER, INTENT (in) :: kmax ! number of levels |
---|
| 945 | INTEGER, INTENT(in) :: ftype ! type of stretching function |
---|
| 946 | REAL(wp) :: pf1 ! "s" stretched value |
---|
| 947 | ! SF12 stretch. funct. parameters |
---|
| 948 | REAL(wp) :: psmth ! smoothing parameter for transition |
---|
| 949 | ! between shallow and deep areas |
---|
| 950 | REAL(wp) :: za1,za2,za3 ! local variables |
---|
| 951 | REAL(wp) :: zn1,zn2,ps ! local variables |
---|
| 952 | REAL(wp) :: za,zb,zx ! local variables |
---|
| 953 | !!---------------------------------------------------------------------- |
---|
| 954 | SELECT CASE (ftype) |
---|
| 955 | |
---|
| 956 | CASE (0) ! M96 stretching function |
---|
| 957 | pf1 = ( TANH( ca * ( s + cb ) ) - TANH( cb * ca ) ) & |
---|
| 958 | & * ( COSH( ca ) & |
---|
| 959 | & + COSH( ca * ( 2.e0 * cb - 1.e0 ) ) ) & |
---|
| 960 | & / ( 2._wp * SINH( ca ) ) |
---|
| 961 | CASE (1) ! SH94 stretching function |
---|
| 962 | IF ( ca == 0 ) THEN ! uniform sigma |
---|
| 963 | pf1 = s |
---|
| 964 | ELSE ! stretched sigma |
---|
| 965 | pf1 = (1._wp - cb) * SINH(ca*s) / SINH(ca) + & |
---|
| 966 | & cb * ( ( TANH(ca*(s + 0.5_wp)) - TANH(0.5_wp*ca) ) / & |
---|
| 967 | & (2._wp*TANH(0.5_wp*ca)) ) |
---|
| 968 | END IF |
---|
| 969 | CASE (2) ! SF12 stretching function |
---|
| 970 | psmth = 1.0_wp ! We consider only the case for efold = 0 |
---|
| 971 | ps = -s |
---|
| 972 | zn1 = 1. / REAL(kmax-1) |
---|
| 973 | zn2 = 1. - zn1 |
---|
| 974 | |
---|
| 975 | za1 = (alpha+2.0_wp)*zn1**(alpha+1.0_wp)-(alpha+1.0_wp)*zn1**(alpha+2.0_wp) |
---|
| 976 | za2 = (alpha+2.0_wp)*zn2**(alpha+1.0_wp)-(alpha+1.0_wp)*zn2**(alpha+2.0_wp) |
---|
| 977 | za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1) |
---|
| 978 | |
---|
| 979 | za = cb - za3*(ca-za1)-za2 |
---|
| 980 | za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp)) ) |
---|
| 981 | zb = (ca - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1) |
---|
| 982 | zx = 1.0_wp-za/2.0_wp-zb |
---|
| 983 | |
---|
| 984 | pf1 = za*(ps*(1.0_wp-ps/2.0_wp))+zb*ps**3.0_wp + & |
---|
| 985 | & zx*( (alpha+2.0_wp)*ps**(alpha+1.0_wp) - & |
---|
| 986 | & (alpha+1.0_wp)*ps**(alpha+2.0_wp) ) |
---|
| 987 | pf1 = -pf1*psmth+ps*(1.0_wp-psmth) |
---|
| 988 | CASE (3) ! stretching function for deeper part of the domain |
---|
| 989 | IF ( ca == 0 ) THEN ! uniform sigma |
---|
| 990 | pf1 = s |
---|
| 991 | ELSE ! stretched sigma |
---|
| 992 | pf1 = (1._wp - cb) * SINH(ca*s) / SINH(ca) |
---|
| 993 | END IF |
---|
| 994 | END SELECT |
---|
| 995 | |
---|
| 996 | END FUNCTION s_coord |
---|
| 997 | |
---|
| 998 | ! ===================================================================================================== |
---|
| 999 | |
---|
| 1000 | FUNCTION D1s_coord( s, ca, cb, alpha, kmax, ftype) RESULT( pf1 ) |
---|
| 1001 | !!---------------------------------------------------------------------- |
---|
| 1002 | !! *** ROUTINE D1s_coord *** |
---|
| 1003 | !! |
---|
| 1004 | !! ** Purpose : provide the 1st derivative of the analytical |
---|
| 1005 | !! stretching function for s-coordinate. |
---|
| 1006 | !! |
---|
| 1007 | !! ** Method : if ftype = 2: Siddorn and Furner 2012 |
---|
| 1008 | !! analytical function is used |
---|
| 1009 | !! if ftype = 1: Song and Haidvogel 1994 |
---|
| 1010 | !! analytical function is used |
---|
| 1011 | !! if ftype = 0: Madec et al. 1996 |
---|
| 1012 | !! analytical function is used |
---|
| 1013 | !! (pag 65 of NEMO Manual) |
---|
| 1014 | !! Reference : Madec, Lott, Delecluse and |
---|
| 1015 | !! Crepon, 1996. JPO, 26, 1393-1408 |
---|
| 1016 | !! |
---|
| 1017 | !! s MUST be NEGATIVE: -1 <= s <= 0 |
---|
| 1018 | !!---------------------------------------------------------------------- |
---|
| 1019 | REAL(wp), INTENT(in ) :: s ! continuous not stretched |
---|
| 1020 | ! "s" coordinate |
---|
| 1021 | REAL(wp), INTENT(in) :: ca ! surface stretch. coeff |
---|
| 1022 | REAL(wp), INTENT(in) :: cb ! bottom stretch. coeff |
---|
| 1023 | REAL(wp), INTENT(in) :: alpha ! alpha stretch. coeff for SF12 |
---|
| 1024 | INTEGER, INTENT (in) :: kmax ! number of levels |
---|
| 1025 | INTEGER, INTENT(in ) :: ftype ! type of stretching function |
---|
| 1026 | REAL(wp) :: pf1 ! "s" stretched value |
---|
| 1027 | ! SF12 stretch. funct. parameters |
---|
| 1028 | REAL(wp) :: psmth ! smoothing parameter for transition |
---|
| 1029 | ! between shallow and deep areas |
---|
| 1030 | REAL(wp) :: za1,za2,za3 ! local variables |
---|
| 1031 | REAL(wp) :: zn1,zn2,ps ! local variables |
---|
| 1032 | REAL(wp) :: za,zb,zx ! local variables |
---|
| 1033 | REAL(wp) :: zt1,zt2,zt3 ! local variables |
---|
| 1034 | !!---------------------------------------------------------------------- |
---|
| 1035 | SELECT CASE (ftype) |
---|
| 1036 | |
---|
| 1037 | CASE (0) ! M96 stretching function |
---|
| 1038 | pf1 = ( ca * ( COSH(ca*(2._wp * cb - 1._wp)) + COSH(ca)) * & |
---|
| 1039 | sech(ca * (s+cb))**2 ) / (2._wp * SINH(ca)) |
---|
| 1040 | CASE (1) ! SH94 stretching function |
---|
| 1041 | IF ( ca == 0 ) then ! uniform sigma |
---|
| 1042 | pf1 = 1._wp / REAL(kmax,wp) |
---|
| 1043 | ELSE ! stretched sigma |
---|
| 1044 | pf1 = (1._wp - cb) * ca * COSH(ca*s) / (SINH(ca) * REAL(kmax,wp)) + & |
---|
| 1045 | cb * ca * & |
---|
| 1046 | (sech((s+0.5_wp)*ca)**2) / (2._wp * TANH(0.5_wp*ca) * REAL(kmax,wp)) |
---|
| 1047 | END IF |
---|
| 1048 | CASE (2) ! SF12 stretching function |
---|
| 1049 | ps = -s |
---|
| 1050 | zn1 = 1. / REAL(kmax-1) |
---|
| 1051 | zn2 = 1. - zn1 |
---|
| 1052 | |
---|
| 1053 | za1 = (alpha+2.0_wp)*zn1**(alpha+1.0_wp)-(alpha+1.0_wp)*zn1**(alpha+2.0_wp) |
---|
| 1054 | za2 = (alpha+2.0_wp)*zn2**(alpha+1.0_wp)-(alpha+1.0_wp)*zn2**(alpha+2.0_wp) |
---|
| 1055 | za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1) |
---|
| 1056 | |
---|
| 1057 | za = cb - za3*(ca-za1)-za2 |
---|
| 1058 | za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp)) ) |
---|
| 1059 | zb = (ca - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1) |
---|
| 1060 | zx = (alpha+2.0_wp)*(alpha+1.0_wp) |
---|
| 1061 | |
---|
| 1062 | zt1 = 0.5_wp*za*(ps-1._wp)*((ps**2.0_wp + 3._wp*ps + 2._wp)*ps**alpha - 2._wp) |
---|
| 1063 | zt2 = zb*(zx*ps**(alpha+1.0_wp) - zx*ps**(alpha) + 3._wp*ps**2._wp) |
---|
| 1064 | zt3 = zx*(ps-1._wp)*ps**(alpha) |
---|
| 1065 | |
---|
| 1066 | pf1 = zt1 + zt2 - zt3 |
---|
| 1067 | |
---|
| 1068 | END SELECT |
---|
| 1069 | |
---|
| 1070 | END FUNCTION D1s_coord |
---|
| 1071 | |
---|
| 1072 | ! ===================================================================================================== |
---|
| 1073 | |
---|
| 1074 | FUNCTION z_mes(dep_up, dep_dw, Cs, s, hc) RESULT( z ) |
---|
| 1075 | !!---------------------------------------------------------------------- |
---|
| 1076 | !! *** ROUTINE z_mes *** |
---|
| 1077 | !! |
---|
| 1078 | !! ** Purpose : provide the analytical trasformation from the |
---|
| 1079 | !! computational space (mes-coordinate) to the |
---|
| 1080 | !! physical space (depth z) |
---|
| 1081 | !! |
---|
| 1082 | !! N.B.: z is downward positive defined as well as envelope surfaces. |
---|
| 1083 | !! Therefore, Cs MUST be positive, meaning 0. <= Cs <= 1. |
---|
| 1084 | !!---------------------------------------------------------------------- |
---|
| 1085 | REAL(wp), INTENT(in ) :: dep_up ! shallower envelope |
---|
| 1086 | REAL(wp), INTENT(in ) :: dep_dw ! deeper envelope |
---|
| 1087 | REAL(wp), INTENT(in ) :: Cs ! stretched s coordinate |
---|
| 1088 | REAL(wp), INTENT(in ) :: s ! not stretched s coordinate |
---|
| 1089 | REAL(wp), INTENT(in ) :: hc ! critical depth |
---|
| 1090 | REAL :: z ! downward positive depth |
---|
| 1091 | !!---------------------------------------------------------------------- |
---|
| 1092 | ! |
---|
| 1093 | |
---|
| 1094 | z = dep_up + hc * s + Cs * (dep_dw - hc - dep_up) |
---|
| 1095 | |
---|
| 1096 | END FUNCTION z_mes |
---|
| 1097 | |
---|
| 1098 | ! ===================================================================================================== |
---|
| 1099 | |
---|
| 1100 | FUNCTION D1z_mes(dep_up, dep_dw, d1Cs, hc, kmax) RESULT( d1z ) |
---|
| 1101 | !!---------------------------------------------------------------------- |
---|
| 1102 | !! *** ROUTINE D1z_mes *** |
---|
| 1103 | !! |
---|
| 1104 | !! ** Purpose : provide the 1st derivative of the analytical |
---|
| 1105 | !! trasformation from the computational space |
---|
| 1106 | !! (mes-coordinate) to the physical space (depth z) |
---|
| 1107 | !! |
---|
| 1108 | !! N.B.: z is downward positive defined as well as envelope surfaces. |
---|
| 1109 | !! Therefore, Cs MUST be positive, meaning 0. <= Cs <= 1. |
---|
| 1110 | !!---------------------------------------------------------------------- |
---|
| 1111 | REAL(wp), INTENT(in ) :: dep_up ! shallower envelope |
---|
| 1112 | REAL(wp), INTENT(in ) :: dep_dw ! deeper envelope |
---|
| 1113 | REAL(wp), INTENT(in ) :: d1Cs ! 1st derivative of the |
---|
| 1114 | ! stretched s coordinate |
---|
| 1115 | REAL(wp), INTENT(in ) :: hc ! critical depth |
---|
| 1116 | INTEGER, INTENT(in ) :: kmax ! number of levels |
---|
| 1117 | REAL(wp) :: d1z ! 1st derivative of downward |
---|
| 1118 | ! positive depth |
---|
| 1119 | !!---------------------------------------------------------------------- |
---|
| 1120 | ! |
---|
| 1121 | |
---|
| 1122 | d1z = hc / REAL(kmax,wp) + d1Cs * (dep_dw - hc - dep_up) |
---|
| 1123 | !d1z = hc + d1Cs * (dep_dw - hc - dep_up) |
---|
| 1124 | |
---|
| 1125 | END FUNCTION D1z_mes |
---|
| 1126 | |
---|
| 1127 | ! ===================================================================================================== |
---|
| 1128 | |
---|
| 1129 | FUNCTION cub_spl(tau, d, n, ibcbeg, ibcend) RESULT ( c ) |
---|
| 1130 | !!---------------------------------------------------------------------- |
---|
| 1131 | !! |
---|
| 1132 | !! CUBSPL defines an interpolatory cubic spline. |
---|
| 1133 | !! |
---|
| 1134 | !! Discussion: |
---|
| 1135 | !! |
---|
| 1136 | !! A tridiagonal linear system for the unknown slopes S(I) of |
---|
| 1137 | !! F at TAU(I), I=1,..., N, is generated and then solved by Gauss |
---|
| 1138 | !! elimination, with S(I) ending up in C(2,I), for all I. |
---|
| 1139 | !! |
---|
| 1140 | !! Author: Carl de Boor |
---|
| 1141 | !! |
---|
| 1142 | !! Reference: Carl de Boor, Practical Guide to Splines, |
---|
| 1143 | !! Springer, 2001, ISBN: 0387953663, LC: QA1.A647.v27. |
---|
| 1144 | !! |
---|
| 1145 | !! Parameters: |
---|
| 1146 | !! |
---|
| 1147 | !! Input: |
---|
| 1148 | !! TAU(N) : the abscissas or X values of the data points. |
---|
| 1149 | !! The entries of TAU are assumed to be strictly |
---|
| 1150 | !! increasing. |
---|
| 1151 | !! |
---|
| 1152 | !! N : the number of data points. N is assumed to be |
---|
| 1153 | !! at least 2. |
---|
| 1154 | !! |
---|
| 1155 | !! IBCBEG : boundary condition indicator at TAU(1). |
---|
| 1156 | !! = 0 no boundary condition at TAU(1) is given. |
---|
| 1157 | !! In this case, the "not-a-knot condition" |
---|
| 1158 | !! is used. |
---|
| 1159 | !! = 1 the 1st derivative at TAU(1) is equal to |
---|
| 1160 | !! the input value D(2,1). |
---|
| 1161 | !! = 2 the 2nd derivative at TAU(1) is equal to |
---|
| 1162 | !! the input value D(2,1). |
---|
| 1163 | !! |
---|
| 1164 | !! IBCEND : boundary condition indicator at TAU(N). |
---|
| 1165 | !! = 0 no boundary condition at TAU(N) is given. |
---|
| 1166 | !! In this case, the "not-a-knot condition" |
---|
| 1167 | !! is used. |
---|
| 1168 | !! = 1 the 1st derivative at TAU(N) is equal to |
---|
| 1169 | !! the input value D(2,2). |
---|
| 1170 | !! = 2 the 2nd derivative at TAU(N) is equal to |
---|
| 1171 | !! the input value D(2,2). |
---|
| 1172 | !! |
---|
| 1173 | !! D(1,1): value of the function at TAU(1) |
---|
| 1174 | !! D(1,1): value of the function at TAU(N) |
---|
| 1175 | !! D(2,1): if IBCBEG is 1 (2) it is the value of the |
---|
| 1176 | !! 1st (2nd) derivative at TAU(1) |
---|
| 1177 | !! D(2,2): if IBCBEG is 1 (2) it is the value of the |
---|
| 1178 | !! 1st (2nd) derivative at TAU(N) |
---|
| 1179 | !! |
---|
| 1180 | !! Output: |
---|
| 1181 | !! C(4,N) : contains the polynomial coefficients of the |
---|
| 1182 | !! cubic interpolating spline. |
---|
| 1183 | !! In the interval interval (TAU(I), TAU(I+1)), |
---|
| 1184 | !! the spline F is given by |
---|
| 1185 | !! |
---|
| 1186 | !! F(X) = |
---|
| 1187 | !! C(1,I) + |
---|
| 1188 | !! C(2,I) * H + |
---|
| 1189 | !! C(3,I) * H^2 / 2 + |
---|
| 1190 | !! C(4,I) * H^3 / 6. |
---|
| 1191 | !! |
---|
| 1192 | !! where H = X - TAU(I). |
---|
| 1193 | !! |
---|
| 1194 | !!---------------------------------------------------------------------- |
---|
| 1195 | INTEGER, INTENT(in ) :: n |
---|
| 1196 | REAL(wp), INTENT(in ) :: tau(n) |
---|
| 1197 | REAL(wp), INTENT(in ) :: d(2,2) |
---|
| 1198 | INTEGER, INTENT(in ) :: ibcbeg, ibcend |
---|
| 1199 | REAL(wp) :: c(4,n) |
---|
| 1200 | REAL(wp) :: divdf1 |
---|
| 1201 | REAL(wp) :: divdf3 |
---|
| 1202 | REAL(wp) :: dtau |
---|
| 1203 | REAL(wp) :: g |
---|
| 1204 | INTEGER(wp) :: i |
---|
| 1205 | !!---------------------------------------------------------------------- |
---|
| 1206 | ! Initialise c and copy d values to c |
---|
| 1207 | c(:,:) = 0.0D+00 |
---|
| 1208 | c(1,1) = d(1,1) |
---|
| 1209 | c(1,n) = d(1,2) |
---|
| 1210 | c(2,1) = d(2,1) |
---|
| 1211 | c(2,n) = d(2,2) |
---|
| 1212 | ! |
---|
| 1213 | ! C(3,*) and C(4,*) are used initially for temporary storage. |
---|
| 1214 | ! Store first differences of the TAU sequence in C(3,*). |
---|
| 1215 | ! Store first divided difference of data in C(4,*). |
---|
| 1216 | DO i = 2, n |
---|
| 1217 | c(3,i) = tau(i) - tau(i-1) |
---|
| 1218 | c(4,i) = ( c(1,i) - c(1,i-1) ) / c(3,i) |
---|
| 1219 | END DO |
---|
| 1220 | |
---|
| 1221 | ! Construct the first equation from the boundary condition |
---|
| 1222 | ! at the left endpoint, of the form: |
---|
| 1223 | ! |
---|
| 1224 | ! C(4,1) * S(1) + C(3,1) * S(2) = C(2,1) |
---|
| 1225 | ! |
---|
| 1226 | ! IBCBEG = 0: Not-a-knot |
---|
| 1227 | IF ( ibcbeg == 0 ) THEN |
---|
| 1228 | |
---|
| 1229 | IF ( n <= 2 ) THEN |
---|
| 1230 | c(4,1) = 1._wp |
---|
| 1231 | c(3,1) = 1._wp |
---|
| 1232 | c(2,1) = 2._wp * c(4,2) |
---|
| 1233 | ELSE |
---|
| 1234 | c(4,1) = c(3,3) |
---|
| 1235 | c(3,1) = c(3,2) + c(3,3) |
---|
| 1236 | c(2,1) = ( ( c(3,2) + 2._wp * c(3,1) ) * c(4,2) & |
---|
| 1237 | * c(3,3) + c(3,2)**2 * c(4,3) ) / c(3,1) |
---|
| 1238 | END IF |
---|
| 1239 | |
---|
| 1240 | ! IBCBEG = 1: derivative specified. |
---|
| 1241 | ELSE IF ( ibcbeg == 1 ) then |
---|
| 1242 | |
---|
| 1243 | c(4,1) = 1._wp |
---|
| 1244 | c(3,1) = 0._wp |
---|
| 1245 | |
---|
| 1246 | ! IBCBEG = 2: Second derivative prescribed at left end. |
---|
| 1247 | ELSE IF ( ibcbeg == 2 ) then |
---|
| 1248 | |
---|
| 1249 | c(4,1) = 2._wp |
---|
| 1250 | c(3,1) = 1._wp |
---|
| 1251 | c(2,1) = 3._wp * c(4,2) - c(3,2) / 2._wp * c(2,1) |
---|
| 1252 | |
---|
| 1253 | ELSE |
---|
| 1254 | WRITE(ctlmes,*) 'CUBSPL - Error, invalid IBCBEG input option!' |
---|
| 1255 | CALL ctl_stop( ctlmes ) |
---|
| 1256 | END IF |
---|
| 1257 | |
---|
| 1258 | ! If there are interior knots, generate the corresponding |
---|
| 1259 | ! equations and carry out the forward pass of Gauss, |
---|
| 1260 | ! elimination after which the I-th equation reads: |
---|
| 1261 | ! |
---|
| 1262 | ! C(4,I) * S(I) + C(3,I) * S(I+1) = C(2,I). |
---|
| 1263 | |
---|
| 1264 | IF ( n > 2 ) THEN |
---|
| 1265 | |
---|
| 1266 | DO i = 2, n-1 |
---|
| 1267 | g = -c(3,i+1) / c(4,i-1) |
---|
| 1268 | c(2,i) = g * c(2,i-1) + 3._wp * & |
---|
| 1269 | ( c(3,i) * c(4,i+1) + c(3,i+1) * c(4,i) ) |
---|
| 1270 | c(4,i) = g * c(3,i-1) + 2._wp * ( c(3,i) + c(3,i+1)) |
---|
| 1271 | END DO |
---|
| 1272 | |
---|
| 1273 | ! Construct the last equation from the second boundary, |
---|
| 1274 | ! condition of the form |
---|
| 1275 | ! |
---|
| 1276 | ! -G * C(4,N-1) * S(N-1) + C(4,N) * S(N) = C(2,N) |
---|
| 1277 | ! |
---|
| 1278 | ! If 1st der. is prescribed at right end ( ibcend == 1 ), |
---|
| 1279 | ! one can go directly to back-substitution, since the C |
---|
| 1280 | ! array happens to be set up just right for it at this |
---|
| 1281 | ! point. |
---|
| 1282 | |
---|
| 1283 | IF ( ibcend < 1 ) THEN |
---|
| 1284 | ! Not-a-knot and 3 <= N, and either 3 < N or also |
---|
| 1285 | ! not-a-knot at left end point. |
---|
| 1286 | IF ( n /= 3 .OR. ibcbeg /= 0 ) THEN |
---|
| 1287 | g = c(3,n-1) + c(3,n) |
---|
| 1288 | c(2,n) = ( ( c(3,n) + 2._wp * g ) * c(4,n) * c(3,n-1) & |
---|
| 1289 | + c(3,n)**2 * ( c(1,n-1) - c(1,n-2) ) / & |
---|
| 1290 | c(3,n-1) ) / g |
---|
| 1291 | g = - g / c(4,n-1) |
---|
| 1292 | c(4,n) = c(3,n-1) |
---|
| 1293 | c(4,n) = c(4,n) + g * c(3,n-1) |
---|
| 1294 | c(2,n) = ( g * c(2,n-1) + c(2,n) ) / c(4,n) |
---|
| 1295 | ELSE |
---|
| 1296 | ! N = 3 and not-a-knot also at left. |
---|
| 1297 | c(2,n) = 2._wp * c(4,n) |
---|
| 1298 | c(4,n) = 1._wp |
---|
| 1299 | g = -1._wp / c(4,n-1) |
---|
| 1300 | c(4,n) = c(4,n) - c(3,n-1) / c(4,n-1) |
---|
| 1301 | c(2,n) = ( g * c(2,n-1) + c(2,n) ) / c(4,n) |
---|
| 1302 | END IF |
---|
| 1303 | |
---|
| 1304 | ELSE IF ( ibcend == 2 ) THEN |
---|
| 1305 | ! IBCEND = 2: Second derivative prescribed at right endpoint. |
---|
| 1306 | c(2,n) = 3._wp * c(4,n) + c(3,n) / 2._wp * c(2,n) |
---|
| 1307 | c(4,n) = 2._wp |
---|
| 1308 | g = -1._wp / c(4,n-1) |
---|
| 1309 | c(4,n) = c(4,n) - c(3,n-1) / c(4,n-1) |
---|
| 1310 | c(2,n) = ( g * c(2,n-1) + c(2,n) ) / c(4,n) |
---|
| 1311 | END IF |
---|
| 1312 | |
---|
| 1313 | ELSE |
---|
| 1314 | ! N = 2 (assumed to be at least equal to 2!). |
---|
| 1315 | |
---|
| 1316 | IF ( ibcend == 2 ) THEN |
---|
| 1317 | |
---|
| 1318 | c(2,n) = 3._wp * c(4,n) + c(3,n) / 2._wp * c(2,n) |
---|
| 1319 | c(4,n) = 2._wp |
---|
| 1320 | g = -1._wp / c(4,n-1) |
---|
| 1321 | c(4,n) = c(4,n) - c(3,n-1) / c(4,n-1) |
---|
| 1322 | c(2,n) = ( g * c(2,n-1) + c(2,n) ) / c(4,n) |
---|
| 1323 | |
---|
| 1324 | ELSE IF ( ibcend == 0 .AND. ibcbeg /= 0 ) THEN |
---|
| 1325 | |
---|
| 1326 | c(2,n) = 2._wp * c(4,n) |
---|
| 1327 | c(4,n) = 1._wp |
---|
| 1328 | g = -1._wp / c(4,n-1) |
---|
| 1329 | c(4,n) = c(4,n) - c(3,n-1) / c(4,n-1) |
---|
| 1330 | c(2,n) = ( g * c(2,n-1) + c(2,n) ) / c(4,n) |
---|
| 1331 | |
---|
| 1332 | ELSE IF ( ibcend == 0 .AND. ibcbeg == 0 ) THEN |
---|
| 1333 | |
---|
| 1334 | c(2,n) = c(4,n) |
---|
| 1335 | |
---|
| 1336 | END IF |
---|
| 1337 | |
---|
| 1338 | END IF |
---|
| 1339 | |
---|
| 1340 | ! Back solve the upper triangular system |
---|
| 1341 | ! |
---|
| 1342 | ! C(4,I) * S(I) + C(3,I) * S(I+1) = B(I) |
---|
| 1343 | ! |
---|
| 1344 | ! for the slopes C(2,I), given that S(N) is already known. |
---|
| 1345 | ! |
---|
| 1346 | DO i = n-1, 1, -1 |
---|
| 1347 | c(2,i) = ( c(2,i) - c(3,i) * c(2,i+1) ) / c(4,i) |
---|
| 1348 | END DO |
---|
| 1349 | |
---|
| 1350 | ! Generate cubic coefficients in each interval, that is, the |
---|
| 1351 | ! derivatives at its left endpoint, from value and slope at its |
---|
| 1352 | ! endpoints. |
---|
| 1353 | DO i = 2, n |
---|
| 1354 | dtau = c(3,i) |
---|
| 1355 | divdf1 = ( c(1,i) - c(1,i-1) ) / dtau |
---|
| 1356 | divdf3 = c(2,i-1) + c(2,i) - 2._wp * divdf1 |
---|
| 1357 | c(3,i-1) = 2._wp * ( divdf1 - c(2,i-1) - divdf3 ) / dtau |
---|
| 1358 | c(4,i-1) = 6._wp * divdf3 / dtau**2 |
---|
| 1359 | END DO |
---|
| 1360 | |
---|
| 1361 | END FUNCTION cub_spl |
---|
| 1362 | |
---|
| 1363 | ! ===================================================================================================== |
---|
| 1364 | |
---|
| 1365 | FUNCTION z_cub_spl(s, tau, c) RESULT ( z_cs ) |
---|
| 1366 | !---------------------------------------------------------------------- |
---|
| 1367 | ! *** function z_cub_spl *** |
---|
| 1368 | ! |
---|
| 1369 | ! ** Purpose : evaluate the cubic spline F in the interval |
---|
| 1370 | ! (TAU(I), TAU(I+1)). F is given by |
---|
| 1371 | ! |
---|
| 1372 | ! |
---|
| 1373 | ! F(X) = C1 + C2*H + (C3*H^2)/2 + (C4*H^3)/6 |
---|
| 1374 | ! |
---|
| 1375 | ! where H = S-TAU(I). |
---|
| 1376 | ! |
---|
| 1377 | ! s MUST be positive: 0 <= s <= 1 |
---|
| 1378 | !--------------------------------------------------------------------- |
---|
| 1379 | INTEGER, PARAMETER :: n = 2 |
---|
| 1380 | REAL(wp), INTENT(in ) :: s |
---|
| 1381 | REAL(wp), INTENT(in ) :: tau(n) |
---|
| 1382 | REAL(wp), INTENT(in ) :: c(4,n) |
---|
| 1383 | REAL(wp) :: z_cs |
---|
| 1384 | REAL(wp) :: c1, c2, c3, c4 |
---|
| 1385 | REAL(wp) :: H |
---|
| 1386 | !--------------------------------------------------------------------- |
---|
| 1387 | c1 = c(1,1) |
---|
| 1388 | c2 = c(2,1) |
---|
| 1389 | c3 = c(3,1) |
---|
| 1390 | c4 = c(4,1) |
---|
| 1391 | |
---|
| 1392 | H = s - tau(1) |
---|
| 1393 | |
---|
| 1394 | z_cs = c1 + c2 * H + (c3 * H**2._wp)/2._wp + (c4 * H**3._wp)/6._wp |
---|
| 1395 | |
---|
| 1396 | END FUNCTION z_cub_spl |
---|
| 1397 | |
---|
| 1398 | ! ===================================================================================================== |
---|
| 1399 | |
---|
| 1400 | FUNCTION D1z_cub_spl(s, tau, c, kmax) RESULT ( d1z_cs ) |
---|
| 1401 | !---------------------------------------------------------------------- |
---|
| 1402 | ! *** function D1z_cub_spl *** |
---|
| 1403 | ! |
---|
| 1404 | ! ** Purpose : evaluate the 1st derivative of cubic spline F |
---|
| 1405 | ! in the interval (TAU(I), TAU(I+1)). |
---|
| 1406 | ! The 1st derivative D1F is given by |
---|
| 1407 | ! |
---|
| 1408 | ! |
---|
| 1409 | ! D1F(S) = C1 + C2*H + (C3*H^2)/2 + (C4*H^3)/6 |
---|
| 1410 | ! |
---|
| 1411 | ! where H = S-TAU(I). |
---|
| 1412 | ! |
---|
| 1413 | ! s MUST be positive: 0 <= s <= 1 |
---|
| 1414 | !--------------------------------------------------------------------- |
---|
| 1415 | INTEGER, PARAMETER :: n = 2 |
---|
| 1416 | REAL(wp), INTENT(in ) :: s |
---|
| 1417 | REAL(wp), INTENT(in ) :: tau(n) |
---|
| 1418 | REAL(wp), INTENT(in ) :: c(4,n) |
---|
| 1419 | INTEGER, INTENT(in ) :: kmax |
---|
| 1420 | REAL(wp) :: d1z_cs |
---|
| 1421 | REAL(wp) :: c2, c3, c4 |
---|
| 1422 | REAL(wp) :: H |
---|
| 1423 | !--------------------------------------------------------------------- |
---|
| 1424 | c2 = c(2,1) / REAL(kmax,wp) |
---|
| 1425 | c3 = c(3,1) / REAL(kmax,wp) |
---|
| 1426 | c4 = c(4,1) / REAL(kmax,wp) |
---|
| 1427 | |
---|
| 1428 | H = s - tau(1) |
---|
| 1429 | |
---|
| 1430 | d1z_cs = c2 + c3 * H + (c4 * H**2._wp)/2._wp |
---|
| 1431 | |
---|
| 1432 | END FUNCTION D1z_cub_spl |
---|
| 1433 | |
---|
| 1434 | ! ===================================================================================================== |
---|
| 1435 | |
---|
| 1436 | |
---|
| 1437 | END MODULE mes |
---|