Changeset 15134
- Timestamp:
- 2021-07-22T18:26:16+02:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/mes.F90
r15131 r15134 4 4 !! Ocean initialization : Multiple Enveloped s coordinate (MES) 5 5 !!============================================================================== 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 7 8 !!---------------------------------------------------------------------- 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 ! 14 10 USE oce ! ocean variables 15 11 USE dom_oce ! ocean domain … … 35 31 INTEGER :: nn_strt(max_nn_env) ! Array specifing the stretching function for each envelope: 36 32 ! 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.)38 33 REAL(wp) :: max_dep_env(max_nn_env) ! global maximum depth of envelopes 39 34 REAL(wp) :: min_dep_env(max_nn_env) ! global minimum depth of envelopes … … 67 62 !! 68 63 !! ** 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: 73 70 !! 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 95 77 !! 96 78 !! Lines represent W-levels: … … 138 120 139 121 ! 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 165 152 !!---------------------------------------------------------------------- 166 153 ! … … 169 156 CALL mes_init 170 157 ! 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-points179 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 185 158 ! Initialize mbathy to the maximum ocean level available 186 187 159 mbathy(:,:) = jpkm1 188 160 189 161 ! Defining scosrf and scobot 190 191 162 scosrf(:,:) = 0._wp ! ocean surface depth (here zero: no under ice-shelf sea) 192 163 scobot(:,:) = bathy(:,:) ! ocean bottom depth 193 164 194 !==================================================================195 ! 1. Defining 1D MEs-levels and computing their depths196 ! (they will be used in diawri.F90)197 !==================================================================198 199 ! Computing gdept_1d and gdepw_1d200 201 DO je = 1, tot_env ! LOOP over all the requested envelopes202 IF (je == 1) THEN203 s_1st = 1204 num_s = nn_slev(je)205 dep_up = 0._wp ! surface206 ELSE207 s_1st = s_1st + num_s - 1208 num_s = nn_slev(je) + 1209 dep_up = rn_ebot_max(je-1)210 ENDIF211 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) THEN219 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 ENDIF223 224 IF ( isodd(je) ) THEN225 DO jk = 1, num_s226 ind = s_1st+jk-1227 ! W-GRID228 z_gsigw1(ind) = -s_coord(sigma(1,jk,num_s), coefa, &229 coefb, alpha, num_s, nn_strt(je))230 ! T-GRID231 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) THEN237 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 ELSE242 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 IF247 ENDDO248 !249 ELSE250 ! Ensure that for splines we do not use SF12251 !IF (nn_strt(je) > 1) nn_strt(je) = 1252 !253 ! Interval's endpoints TAU are 0 and 1.254 tau(1) = 0._wp255 tau(n) = 1._wp256 257 IF ( je == 2 ) THEN258 dep0 = 0._wp259 ELSE260 dep0 = rn_ebot_max(je-2)261 END IF262 dep1 = dep_up263 dep2 = dep_dw264 IF ( je < tot_env ) dep3 = rn_ebot_max(je+1)265 !266 ! Computing derivative analytically267 !268 ! 1. Derivative for upper envelope269 ibcbeg = 1270 IF (nn_strt(je-1) == 2) THEN271 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 ELSE277 alpha = 0._wp278 num_s_up = 1279 coefa = rn_e_th(je-1)280 coefb = rn_e_bb(je-1)281 END IF282 D1s_up = D1s_coord(-1._wp, coefa, coefb, &283 alpha, num_s_up, nn_strt(je-1))284 d(1,1) = dep_up285 IF (nn_strt(je-1) == 2) THEN286 d(2,1) = D1z_mes(dep0, dep1, D1s_up, 0._wp, 1)287 ELSE288 d(2,1) = D1z_mes(dep0, dep1, D1s_up, rn_e_hc(je-1), 1)289 END IF290 !291 ! 2. Derivative for deeper envelope292 IF ( je < tot_env ) THEN293 ibcend = 1 ! Bound. cond for the 1st derivative294 IF (nn_strt(je+1) == 2) THEN295 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 ELSE301 alpha = 0._wp302 num_s_dw = 1303 coefa = rn_e_th(je+1)304 coefb = rn_e_bb(je+1)305 END IF306 D1s_dw = D1s_coord(0._wp, coefa, coefb, &307 alpha, num_s_dw, nn_strt(je+1))308 d(1,2) = dep_dw309 IF (nn_strt(je+1) == 2) THEN310 d(2,2) = D1z_mes(dep2, dep3, D1s_dw, 0._wp, 1)311 ELSE312 d(2,2) = D1z_mes(dep2, dep3, D1s_dw, rn_e_hc(je+1), 1)313 END IF314 ELSE315 ibcend = 2 ! Bound. cond for the 2nd derivative316 d(2,2) = 0._wp ! 2nd der = 0317 ENDIF318 !319 ! Computing spline's coefficients320 IF ( lwp ) THEN321 WRITE(numout,*) "BOUNDARY COND. to CONSTRAIN SPLINES:"322 WRITE(numout,*) "ibcbeg: ", ibcbeg323 WRITE(numout,*) "ibcend: ", ibcend324 ENDIF325 c = cub_spl(tau, d, n, ibcbeg, ibcend )326 327 IF ( lwp ) THEN328 WRITE(numout,*) "Subzone between depth ", dep1, " and depth ", dep2329 WRITE(numout,*) ""330 WRITE(numout,*) "D1s_up = ", D1s_up331 WRITE(numout,*) "1st deriv. of z(mes) at depth ", dep1, " is ", d(2,1)332 WRITE(numout,*) "D1s_dw = ", D1s_dw333 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 IF350 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_s357 ind = s_1st+jk-1358 ! W-GRID359 z_gsigw1(s_1st+jk-1) = -s_coord( sigma(1,jk,num_s), coefa, &360 coefb, alpha, num_s, nn_strt(je) )361 ! T-GRID362 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 DO369 370 END IF371 372 END DO373 374 !==================================================================375 ! 2. Computing derivative of MEs-coordinate and scale factors376 ! as finite differences377 !==================================================================378 379 ! Computing derivative of MEs-coordinate as finite differences380 DO jk = 1, jpkm1381 z_esigt1(jk) = z_gsigw1(jk+1) - z_gsigw1(jk)382 z_esigw1(jk+1) = z_gsigt1(jk+1) - z_gsigt1(jk)383 END DO384 ! Surface385 z_esigw1( 1 ) = 2._wp * ( z_gsigt1( 1 ) - z_gsigw1( 1 ) )386 ! Bottom387 z_esigt1(jpk) = 2._wp * ( z_gsigt1(jpk) - z_gsigw1(jpk) )388 389 ! Computing scale factors as finite differences390 DO jk=1,jpkm1391 e3t_1d(jk) = gdepw_1d(jk+1) - gdepw_1d(jk)392 e3w_1d(jk+1) = gdept_1d(jk+1) - gdept_1d(jk)393 ENDDO394 ! Bottom:395 jk = jpk396 e3t_1d(jk) = 2.0_wp * (gdept_1d(jk) - gdepw_1d(jk))397 ! and Surface (H. Liu)398 jk = 1399 e3w_1d(jk) = 2.0_wp * (gdept_1d(1) - gdepw_1d(1))400 401 ! Control print402 IF ( lwp ) THEN403 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, jpk409 WRITE(numout,*) jk, z_gsigw1(jk), z_esigw1(jk), z_gsigt1(jk), z_esigt1(jk)410 ENDDO411 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, jpk419 WRITE(numout,*) jk, gdepw_1d(jk), e3w_1d(jk), gdept_1d(jk), e3t_1d(jk)420 ENDDO421 WRITE(numout,*) "-----------------------------------------------------------------"422 ENDIF423 424 ! Checking monotonicity for cubic splines425 DO jk = 1, jpk-1426 IF ( gdept_1d(jk+1) < gdept_1d(jk) ) THEN427 WRITE(ctlmes,*) 'NOT MONOTONIC gdept_1d: change envelopes'428 CALL ctl_stop( ctlmes )429 END IF430 IF ( gdepw_1d(jk+1) < gdepw_1d(jk) ) THEN431 WRITE(ctlmes,*) 'NOT MONOTONIC gdepw_1d: change envelopes'432 CALL ctl_stop( ctlmes )433 END IF434 END DO435 436 165 !======================== 437 ! 3.Compute 3D MEs mesh166 ! Compute 3D MEs mesh 438 167 !======================== 439 168 … … 449 178 env_up = envlt(:,:,je-1) 450 179 ENDIF 451 452 180 env_dw = envlt(:,:,je) 453 env_hc = rn_e_hc(je)454 alpha = rn_e_al(je)455 181 456 182 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) 457 188 ! 458 189 DO jk = 1, num_s 459 190 ind = s_1st+jk-1 191 ! 460 192 DO jj = 1,jpj 461 193 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 465 203 zcoeft = -sigma(0, jk, num_s) 466 204 !IF (nn_strt(je) == 2) THEN … … 471 209 ! zcoeft = zcoeft*(env_hc/(env_dw(ji,jj)-env_up(ji,jj))) 472 210 !ENDIF 473 ELSE ! deep water, stretched s 211 ! deep water, stretched s 212 ELSE 474 213 IF (nn_strt(je) == 2) THEN 475 214 coefa = rn_e_th(je) / (env_dw(ji,jj)-env_up(ji,jj)) 476 215 coefb = rn_e_bb(je) + ((env_dw(ji,jj)-env_up(ji,jj))*rn_e_ba(je)) 477 216 coefb = 1.0_wp-(coefb/(env_dw(ji,jj)-env_up(ji,jj))) 478 ! W-GRID479 zcoefw = -s_coord( sigma(1,jk,num_s), coefa, &480 coefb, alpha, num_s, nn_strt(je) )481 ! T-GRID482 zcoeft = -s_coord( sigma(0,jk,num_s), coefa, &483 coefb, alpha, num_s, nn_strt(je) )484 ELSE485 zcoefw = z_gsigw1(ind)486 zcoeft = z_gsigt1(ind)487 217 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) ) 488 224 ENDIF 489 225 ! 226 z_gsigw3(ji,jj,ind) = zcoefw 227 z_gsigt3(ji,jj,ind) = zcoeft 228 ! 229 ! -------------------------------------------------- 230 ! Computing model levels depths 490 231 IF (nn_strt(je) /= 2) THEN 491 232 gdept_0(ji,jj,ind) = z_mes(env_up(ji,jj), env_dw(ji,jj), zcoeft, & … … 591 332 ! WRITE(numout,*) "ibcend: ", ibcend 592 333 !ENDIF 334 593 335 c = cub_spl(tau, d, n, ibcbeg, ibcend ) 594 336 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 596 342 DO jk = 1, num_s 597 343 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) 602 358 END DO 603 359 604 ! Control print605 IF ( lwp ) THEN606 IF ( jj == 25 .AND. ji == 7 ) THEN607 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 IF621 END IF622 623 360 END DO 624 361 END DO … … 644 381 jk = jpk 645 382 e3t_0(:,:,jk) = 2.0_wp * (gdept_0(:,:,jk) - gdepw_0(:,:,jk)) 646 647 IF ( lwp ) THEN648 jj = 25649 ji = 7650 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, jpk654 WRITE(numout,*) jk, gdepw_0(ji,jj,jk), gdept_0(ji,jj,jk), e3w_0(ji,jj,jk), e3t_0(ji,jj,jk)655 ENDDO656 END IF657 658 !==============================659 ! Computing gdept3w660 !==============================661 662 gde3w_0(:,:,1) = 0.5 * e3w_0(:,:,1)663 DO jk = 2, jpk664 gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk)665 END DO666 e3t_0(:,:,jpk) = 2.0_wp * (gdept_0(:,:,jpk) - gdepw_0(:,:,jpk))667 383 668 384 !============================== … … 825 541 !================================================================================ 826 542 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 827 626 DO ji = 1, jpi 828 627 DO jj = 1, jpj 829 DO jk = 1, mbathy(ji,jj)628 DO jk = 1, jpk !mbathy(ji,jj) 830 629 ! check coordinate is monotonically increasing 831 630 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, jk833 WRITE(numout,*) 'ERROR zgr_mes: e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk631 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 834 633 WRITE(numout,*) 'e3w',e3w_0(ji,jj,:) 835 634 WRITE(numout,*) 'e3t',e3t_0(ji,jj,:) … … 838 637 ! and check it has never gone negative 839 638 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, jk841 WRITE(numout,*) 'ERROR zgr_mes: gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk639 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 842 641 WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:) 843 642 WRITE(numout,*) 'gdept',gdept_0(ji,jj,:) 844 643 CALL ctl_stop( ctmp1 ) 845 644 ENDIF 846 ! ! and check it never exceeds the total depth847 ! IF ( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN848 ! WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk849 ! WRITE(numout,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk850 ! WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:)851 ! CALL ctl_stop( ctmp1 )852 ! ENDIF853 645 END DO 854 855 ! DO jk = 1, mbathy(ji,jj)-1856 ! ! and check it never exceeds the total depth857 ! IF ( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN858 ! WRITE(ctmp1,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk859 ! WRITE(numout,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk860 ! WRITE(numout,*) 'gdept',fsdept(ji,jj,:)861 ! CALL ctl_stop( ctmp1 )862 ! ENDIF863 ! END DO864 646 END DO 865 647 END DO 866 648 ! 867 649 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 )870 650 ! 871 651 IF( nn_timing == 1 ) CALL timing_stop('mes_build') … … 876 656 877 657 SUBROUTINE mes_init 658 !!----------------------------------------------------------------------------- 659 !! *** ROUTINE mes_init *** 660 !! 661 !! ** Purpose : Initilaise a Multi Enveloped s-coordinate (MEs) system 662 !! 663 !!----------------------------------------------------------------------------- 664 878 665 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) 885 670 886 671 NAMELIST/namzgr_mes/rn_bot_min , rn_bot_max , ln_envl, & … … 893 678 ! 894 679 CALL wrk_alloc( jpi, jpj, max_nn_env , envlt ) 895 pmsk(:,:) = 0.0896 680 ! 897 681 ! Initialising some arrays 898 rn_ebot_max(:) = 0.0899 682 max_dep_env(:) = 0.0 900 683 min_dep_env(:) = 0.0 … … 963 746 IF (MAXVAL(envlt(:,:,je+1)) < MAXVAL(envlt(:,:,je))) CALL ctl_stop( ctlmes ) 964 747 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 983 749 DO je = 1, tot_env 984 750 max_dep_env(je) = MAXVAL(envlt(:,:,je))
Note: See TracChangeset
for help on using the changeset viewer.
