Changeset 15126
- Timestamp:
- 2021-07-16T17:42:38+02:00 (3 years ago)
- Location:
- NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src
- Files:
-
- 1 added
- 1 edited
-
mes.F90 (added)
-
zgrmes.F90 (modified) (9 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/tools_r4.0-HEAD_dev_MEs/DOMAINcfg/src/zgrmes.F90
r15125 r15126 4 4 !! Ocean initialization : Multiple Enveloped s coordinate (MES) 5 5 !!============================================================================== 6 !! NEMO 3.6 ! 2017-05(D. Bruciaferri)6 !! NEMO 4.0 ! 2021-07 (D. Bruciaferri) 7 7 !!---------------------------------------------------------------------- 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 !!--------------------------------------------------------------------- 8 ! 14 9 USE oce ! ocean variables 15 10 USE dom_oce ! ocean domain 16 11 USE depth_e3 ! depth <=> e3 17 12 USE closea ! closed seas 13 USE mes ! MEs-coordinates 18 14 ! 19 15 USE in_out_manager ! I/O manager … … 28 24 29 25 PUBLIC zgr_mes ! called by domzgr.F90 30 ! Bruciaferri, Shapiro & Wobus (2018) envelopes and stretching parameters 31 CHARACTER(lc) :: ctlmes ! control message error (lc=256) 32 INTEGER, PARAMETER :: max_nn_env = 5 ! Maximum allowed number of envelopes. 33 REAL(wp) :: rn_bot_min ! minimum depth of ocean bottom (>0) (m) 34 REAL(wp) :: rn_bot_max ! maximum depth of ocean bottom (= ocean depth) (>0) (m) 35 LOGICAL :: ln_envl(max_nn_env) ! Array with flags (T/F) specifing which envelope is used 36 INTEGER :: nn_strt(max_nn_env) ! Array specifing the stretching function for each envelope: 37 ! Madec 1996 (0), Song and Haidvogel 1994 (1) 38 REAL(wp) :: rn_ebot_min(max_nn_env) ! minimum depth of envelopes (>0) (m) 39 REAL(wp) :: rn_ebot_max(max_nn_env) ! maximum depth of envelopes (= envelopes depths) (>0) (m) 40 INTEGER :: nn_slev(max_nn_env) ! Array specifing number of levels of each enveloped vertical zone 41 REAL(wp) :: rn_e_hc(max_nn_env) ! Array specifing critical depth for transition to stretched 42 ! coordinates of each envelope 43 REAL(wp) :: rn_e_th(max_nn_env) ! Array specifing surface control parameter (0<=th<=20) of SH94 44 ! or rn_zs of SF12 for each vertical sub-zone 45 REAL(wp) :: rn_e_bb(max_nn_env) ! Array specifing bottom control parameter (0<=bb<=1) of SH94 46 ! or rn_zb_b of SF12 for each vertical sub-zone 47 REAL(wp) :: rn_e_ba(max_nn_env) ! Array specifing bottom control parameter rn_zb_a of SF12 for 48 ! each vertical sub-zone 49 REAL(wp) :: rn_e_al(max_nn_env) ! Array specifing stretching parameter rn_alpha of SF12 for 50 ! each vertical sub-zone 51 LOGICAL :: ln_loc_mes ! To use localised MEs (.TRUE.) or not (.FALSE. 26 ! PROPERTIES OF THE INPUT VERTICAL COORDINATE SYSTEM (ln_loc_mes=.TRUE.) 27 ! 28 REAL(wp), POINTER, DIMENSION(:) :: e3t_1d_in , e3w_1d_in !: ref. scale factors (m) 29 REAL(wp), POINTER, DIMENSION(:) :: gdept_1d_in, gdepw_1d_in !: ref. depth (m) 30 ! 31 REAL(wp), POINTER, DIMENSION(:,:) :: mbathy_in 32 ! 33 REAL(wp), POINTER, DIMENSION(:,:,:) :: e3t_in, e3u_in , e3v_in , e3f_in !: scale factors [m] 34 REAL(wp), POINTER, DIMENSION(:,:,:) :: e3w_in, e3uw_in, e3vw_in !: - - 35 ! 36 REAL(wp), POINTER, DIMENSION(:,:,:) :: gdept_in, gdepw_in !: depths [m] 52 37 53 38 !! * Substitutions … … 68 53 ! PROPERTIES OF THE INPUT VERTICAL COORDINATE SYSTEM (ln_loc_mes=.TRUE.) 69 54 ! 70 REAL(wp), POINTER, DIMENSION(:) :: e3t_1d_in , e3w_1d_in !: ref. scale factors (m)71 REAL(wp), POINTER, DIMENSION(:) :: gdept_1d_in, gdepw_1d_in !: ref. depth (m)72 !73 REAL(wp), POINTER, DIMENSION(:,:,:) :: e3t_in, e3u_in , e3v_in , e3f_in !: scale factors [m]74 REAL(wp), POINTER, DIMENSION(:,:,:) :: e3w_in, e3uw_in, e3vw_in !: - -75 !76 REAL(wp), POINTER, DIMENSION(:,:,:) :: gdept_in, gdepw_in !: depths [m]77 55 !----------------------------------------------------------------------- 78 56 ! … … 80 58 CALL zgr_mes_build 81 59 82 if( ln_loc_mes ) THEN 83 ! Local MEs 60 ! Local MEs 61 IF( ln_loc_mes ) THEN 62 ! 84 63 CALL wrk_alloc( jpk, gdept_1d_in, gdepw_1d_in, e3t_1d_in , e3w_1d_in) 64 CALL wrk_alloc( jpi, jpj, mbathy_in) 85 65 CALL wrk_alloc( jpi, jpj, jpk, gdept_in, gdepw_in, e3t_in, e3w_in) 86 66 CALL wrk_alloc( jpi, jpj, jpk, e3u_in , e3v_in , e3f_in, e3uw_in, e3vw_in) 87 67 ! 88 68 ! Reading the input vertical grid that will be used globally 89 CALL zgr_dom_read(ln_isfcav, & ! ice-cavities 90 & gdept_1d_in, gdepw_1d_in, e3t_1d_in, e3w_1d_in, & ! 1D ref. vert. coord. 91 & gdept_in , gdepw_in , & ! 3D t & w-points depth 92 & e3t_in , e3u_in , e3v_in , e3f_in , & ! vertical scale factors 93 & e3w_in , e3uw_in , e3vw_in ) 69 CALL zgr_dom_read 70 ! Creating the local MEs within the global input vertical grid 71 CALL zgr_mes_local 72 ! 73 CALL wrk_dealloc( jpk, gdept_1d_in, gdepw_1d_in, e3t_1d_in , e3w_1d_in) 74 CALL wrk_dealloc( jpi, jpj, mbathy_in) 75 CALL wrk_dealloc( jpi, jpj, jpk, gdept_in, gdepw_in, e3t_in, e3w_in) 76 CALL wrk_dealloc( jpi, jpj, jpk, e3u_in , e3v_in , e3f_in, e3uw_in, e3vw_in) 77 ! 94 78 END IF 95 79 … … 98 82 ! ===================================================================================================== 99 83 100 SUBROUTINE zgr_mes_build 101 !!----------------------------------------------------------------------------- 102 !! *** ROUTINE zgr_mes *** 103 !! 104 !! ** Purpose : define the Multi Enveloped S-coordinate (MES) system 105 !! 106 !! ** Method : s-coordinate defined respect to multiple 107 !! externally defined envelope 108 !! 109 !! The depth of model levels is defined as the product of an 110 !! analytical function by the externally defined envelopes. 111 !! 112 !! - Read envelopes (in meters) at t-point and compute the 113 !! envelopes at u-, v-, and f-points. 114 !! hbatu = mi( hbatt ) 115 !! hbatv = mj( hbatt ) 116 !! hbatf = mi( mj( hbatt ) ) 117 !! 118 !! - Compute z_gsigt, z_gsigw, z_esigt, z_esigw from an analytical 119 !! function and its finite difference derivative. 120 !! z_gsigt(k) = fssig (k ) 121 !! z_gsigw(k) = fssig (k-0.5) 122 !! z_esigt(k) = z_gsigw(k+1) - z_gsigw(k) 123 !! z_esigw(k) = z_gsigt(k+1) - z_gsigt(k) 124 !! 125 !! Three options for stretching are given: 126 !! 127 !! *) nn_strt = 0 : Madec et al 1996 cosh/tanh function 128 !! *) nn_strt = 1 : Song and Haidvogel 1994 sinh/tanh function 129 !! *) nn_strt = 2 : Siddorn and Furner gamma function 130 !! 131 !!---------------------------------------------------------------------- 132 !! SKETCH of the GEOMETRY OF THE S-S-Z COORDINATE SYSTEM 133 !! 134 !! Lines represent W-levels: 135 !! 136 !! 0: First s-level part. The total number 137 !! of levels in this part is controlled 138 !! by the nn_slev(1) namelist parameter. 139 !! 140 !! Depth first W-lev: 0 m (surfcace) 141 !! Depth last W-lev: depth of 1st envelope 142 !! 143 !! o: Second s-level part. The total number 144 !! of levels in this part is controlled 145 !! by the nn_slev(2) namelist parameter. 146 !! 147 !! Depth last W-lev: depth of 2nd envelope 148 !! 149 !! @: Third s-level part. The total number 150 !! of levels in this part is controlled 151 !! by the nn_slev(3) namelist parameter. 152 !! 153 !! Depth last W-lev: depth of 3rd envelope 154 !! 155 !! z |~~~~~~~~~~~~~~~~~~~~~~~0~~~~~~~~~~~~~~~~~~~~~~~ SURFACE 156 !! | 157 !! |-----------------------0----------------------- ln_envl(1) = true, nn_slev(1) = 3 158 !! | 159 !! |=======================0======================= ENVELOPE 1 160 !! | 161 !! |-----------------------o----------------------- 162 !! | ln_envl(2) = true, nn_slev(2) = 3 163 !! |-----------------------o----------------------- 164 !! | 165 !! |¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬o¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬¬ ENVELOPE 2 166 !! | 167 !! |-----------------------@----------------------- ln_envl(3) = true, nn_slev(3) = 2 168 !! | 169 !! |_______________________@_______________________ ENVELOPE 3 170 !! \|/ 171 !! 172 !! 173 !!----------------------------------------------------------------------------- 174 !! Author: Diego Bruciaferri (University of Plymouth), September 2017 175 !!----------------------------------------------------------------------------- 176 177 ! 178 CHARACTER(lc) :: env_name ! name of the externally defined envelope 179 INTEGER :: ji, jj, jk, jl, je ! dummy loop argument 180 INTEGER :: ios ! Local integer output status for namelist read 181 INTEGER :: inum ! temporary logical unit 182 INTEGER :: tot_env, cor_env 183 INTEGER :: num_s, s_1st, ind ! for loops over envelopes 184 INTEGER :: num_s_up, num_s_dw ! for loops over envelopes 185 REAL(wp) :: D1s_up, D1s_dw ! for loops over envelopes 186 INTEGER :: ibcbeg, ibcend ! boundary condition for cubic splines 187 INTEGER, PARAMETER :: n = 2 ! number of considered points for each interval 188 REAL(wp) :: c(4,n) ! matrix for cubic splines coefficents 189 REAL(wp) :: d(2,2) ! matrix for depth values and depth derivative 190 ! at intervals' boundaries 191 REAL(wp) :: tau(n) ! abscissas of intervals' boundaries 192 REAL(wp) :: coefa, coefb ! for loops over envelopes 193 REAL(wp) :: alpha, env_hc ! for loops over envelopes 194 REAL(wp) :: dep_up, dep_dw ! for loops over envelopes 195 REAL(wp) :: dep0, dep1, dep2, dep3 ! for loops over envelopes 196 ! for cubic splines 197 REAL(wp) :: zcoeft, zcoefw, sw, st ! temporary scalars 198 REAL(wp) :: max_env_up, min_env_up ! to identify geopotential levels 199 REAL(wp) :: max_env_dw, min_env_dw ! to identify geopotential levels 200 REAL(wp), POINTER, DIMENSION(:) :: z_gsigw1, z_gsigt1 ! MES-coordinates analytical trans. 201 ! (0.<= s <= 1.) for T- and W-points 202 REAL(wp), POINTER, DIMENSION(:) :: z_esigw1, z_esigt1 ! MES analytical scale factors 203 ! of T- and W-points 204 REAL(wp), POINTER, DIMENSION(:,: ) :: env_up, env_dw ! for loops over envelopes 205 REAL(wp), POINTER, DIMENSION(:,: ) :: env0, env1, env2, env3 ! for loops over envelopes 206 ! for cubic splines 207 REAL(wp), POINTER, DIMENSION(:,:,:) :: envlt ! array for the envelopes 208 INTEGER :: gst_envl(max_nn_env) ! Array to deal with a ghost last envelope 209 210 NAMELIST/namzgr_mes/rn_bot_min , rn_bot_max , ln_envl, nn_strt, & 211 rn_ebot_min, rn_ebot_max, nn_slev, rn_e_hc, & 212 rn_e_th, rn_e_bb, rn_e_ba, rn_e_al, ln_loc_mes 213 214 !!---------------------------------------------------------------------- 215 ! 216 IF( nn_timing == 1 ) CALL timing_start('zgr_mes') 217 ! 218 CALL wrk_alloc( jpi, jpj, max_nn_env , envlt ) 219 CALL wrk_alloc( jpi, jpj, env_up, env_dw, env0, env1, env2, env3 ) 220 CALL wrk_alloc( jpk, z_gsigw1, z_gsigt1, z_esigw1, z_esigt1 ) 221 ! 222 ! Namelist namzgr_mes in reference namelist : envelopes and 223 ! sigma-stretching parameters 224 REWIND( numnam_ref ) 225 READ ( numnam_ref, namzgr_mes, IOSTAT = ios, ERR = 901) 226 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_mes in reference namelist', lwp ) 227 ! Namelist namzgr_mes in configuration namelist : envelopes and 228 ! sigma-stretching parameters 229 REWIND( numnam_cfg ) 230 READ ( numnam_cfg, namzgr_mes, IOSTAT = ios, ERR = 902 ) 231 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzgr_mes in configuration namelist', lwp ) 232 IF(lwm) WRITE ( numond, namzgr_mes ) 233 234 IF(lwp) THEN ! control print 235 WRITE(numout,*) 236 WRITE(numout,*) 'domzgr_mes : Multi Enveloped S-coordinate (Bruciaferri, Shapiro and Wobus 2017)' 237 WRITE(numout,*) '~~~~~~~~~~~' 238 WRITE(numout,*) ' Namelist namzgr_mes' 239 WRITE(numout,*) '' 240 WRITE(numout,*) ' Minimum depth of the ocean rn_bot_min, ', rn_bot_min 241 WRITE(numout,*) ' Maximum depth of the ocean rn_bot_max, ', rn_bot_max 242 WRITE(numout,*) '' 243 WRITE(numout,*) '-------------------------------------------------------------------------------' 244 DO je = 1, max_nn_env 245 WRITE(numout,*) 'SUBDOMAIN ', je,':' 246 IF ( je == 1) THEN 247 WRITE(numout,*) ' Envelope up : envelope 0, free surface' 248 WRITE(numout,*) ' Envelope down: envelope ',je,',ln_envl(',je,') = ',ln_envl(je) 249 ELSE 250 WRITE(numout,*) ' Envelope up : envelope ',je-1,',ln_envl(',je,') = ',ln_envl(je-1) 251 WRITE(numout,*) ' Envelope down: envelope ',je ,',ln_envl(',je,') = ',ln_envl(je) 252 END IF 253 WRITE(numout,*) ' min dep of envlp down rn_ebot_min(',je,') = ',rn_ebot_min(je) 254 WRITE(numout,*) ' max dep of envlp down rn_ebot_max(',je,') =',rn_ebot_max(je) 255 WRITE(numout,*) ' num. of MEs-lev. nn_slev(',je,') = ',nn_slev(je) 256 IF ( isodd(je) ) THEN 257 WRITE(numout,*) ' Stretched s-coordinates: ' 258 ELSE 259 WRITE(numout,*) ' Stretched CUBIC SPLINES: ' 260 END IF 261 IF (nn_strt(je) == 0) WRITE(numout,*) ' M96 stretching function' 262 IF (nn_strt(je) == 1) WRITE(numout,*) ' SH94 stretching function' 263 IF (nn_strt(je) == 2) WRITE(numout,*) ' SF12 stretching function' 264 WRITE(numout,*) ' critical depth rn_e_hc(',je,') = ',rn_e_hc(je) 265 WRITE(numout,*) ' surface stretc. coef. rn_e_th(',je,') = ',rn_e_th(je) 266 IF (nn_strt(je) == 2) THEN 267 WRITE(numout,*) ' bottom stretc. coef. rn_e_ba(',je,') = ',rn_e_ba(je) 268 END IF 269 WRITE(numout,*) ' bottom stretc. coef. rn_e_bb(',je,') = ',rn_e_bb(je) 270 IF (nn_strt(je) == 2) THEN 271 WRITE(numout,*) ' bottom stretc. coef. rn_e_al(',je,') = ',rn_e_al(je) 272 END IF 273 WRITE(numout,*) '-------------------------------------------------------------------------------' 274 ENDDO 275 ENDIF 276 277 ! Check if namelist is defined correctly. 278 ! Not strictly needed but we force the user 279 ! to define the namelist correctly. 280 281 tot_env = 0 ! total number of requested envelopes 282 cor_env = 0 ! total number of correctly defined envelopes 283 DO je = 1, max_nn_env 284 IF ( ln_envl(je) ) tot_env = tot_env + 1 285 IF ( ln_envl(je) .AND. cor_env == (je-1)) cor_env = cor_env + 1 286 ENDDO 287 WRITE(ctlmes,*) 'number of REQUESTED envelopes and number of CORRECTLY defined envelopes are DIFFERENT' 288 IF ( tot_env /= cor_env ) CALL ctl_stop( ctlmes ) 289 290 !Checking consistency of user defined parameters 291 WRITE(ctlmes,*) 'TOT number of levels (jpk) IS DIFFERENT from sum over nn_slev(:)' 292 IF ( SUM(nn_slev(:)) /= jpk ) CALL ctl_stop( ctlmes ) 293 294 ! Determining if there is a "ghost" envelope: 295 gst_envl(:) = 0 296 IF ( nn_slev(tot_env) == 0 ) gst_envl(tot_env) = 1 297 298 ! Reading bathymetry and envelopes. 299 ! In future should be included in zgr_bat. 300 301 IF( ntopo == 1 ) THEN 302 303 CALL iom_open ( 'bathy_meter.nc', inum ) 304 CALL iom_get ( inum, jpdom_data, 'Bathymetry' , bathy, lrowattr=ln_use_jattr ) 305 306 DO je = 1, tot_env 307 WRITE (env_name, '(A6,I0)') 'hbatt_', je 308 CALL iom_get ( inum, jpdom_data, TRIM(env_name) , envlt(:,:,je), lrowattr=ln_use_jattr ) 309 ENDDO 310 311 CALL iom_close( inum ) 312 313 ELSE 314 315 WRITE(ctlmes,*) 'parameter , ntopo = ', ntopo 316 CALL ctl_stop( ctlmes ) 317 318 ENDIF 319 IF( nn_closea == 0 ) CALL clo_bat( bathy, mbathy ) !== NO closed seas or lakes ==! 320 321 ! Checking consistency of envelopes 322 323 DO je = 1, tot_env-1 324 WRITE(ctlmes,*) 'Envelope ', je+1, ' is shallower that Envelope ', je 325 IF (MAXVAL(envlt(:,:,je+1)) < MAXVAL(envlt(:,:,je))) CALL ctl_stop( ctlmes ) 326 ENDDO 327 328 ! Set maximum and minimum ocean depth 329 bathy(:,:) = MIN( rn_bot_max, bathy(:,:) ) 330 331 DO jj = 1, jpj 332 DO ji = 1, jpi 333 IF( bathy(ji,jj) > 0._wp ) bathy(ji,jj) = MAX( rn_bot_min, bathy(ji,jj) ) 334 END DO 335 END DO 336 337 ! Initializing to 0.0 some arrays: 338 ! 339 ! *) gdept_1d, gdepw_1d : depth of T- and W-point (m) 340 ! *) e3t_1d , e3w_1d : scale factors at T- and W-levels (m) 341 ! *) z_gsigt1, z_gsigw1 : MES-coordinates (0.<= s <= 1.) of T- and W-points 342 343 gdept_1d = 0._wp ; gdepw_1d = 0._wp ; 344 e3t_1d = 0._wp ; e3w_1d = 0._wp ; 345 z_gsigw1 = 0._wp ; z_gsigt1 = 0._wp ; 346 z_esigw1 = 0._wp ; z_esigt1 = 0._wp ; 347 348 ! Initialize mbathy to the maximum ocean level available 349 350 mbathy(:,:) = jpkm1 351 352 ! Defining scosrf and scobot 353 354 scosrf(:,:) = 0._wp ! ocean surface depth (here zero: no under ice-shelf sea) 355 scobot(:,:) = bathy(:,:) ! ocean bottom depth 356 357 !================================================================== 358 ! 1. Defining 1D MEs-levels and computing their depths 359 ! (they will be used in diawri.F90) 360 !================================================================== 361 362 ! Computing gdept_1d and gdepw_1d 363 364 DO je = 1, tot_env ! LOOP over all the requested envelopes 365 IF (je == 1) THEN 366 s_1st = 1 367 num_s = nn_slev(je) 368 dep_up = 0._wp ! surface 369 ELSE 370 s_1st = s_1st + num_s - 1 371 num_s = nn_slev(je) + 1 372 dep_up = rn_ebot_max(je-1) 373 ENDIF 374 375 dep_dw = rn_ebot_max(je) 376 env_hc = rn_e_hc(je) 377 coefa = rn_e_th(je) 378 coefb = rn_e_bb(je) 379 alpha = rn_e_al(je) 380 381 IF (nn_strt(je) == 2) THEN 382 coefa = coefa / (dep_dw - dep_up) 383 coefb = coefb + ((dep_dw - dep_up)*rn_e_ba(je)) 384 coefb = 1.0_wp - (coefb/(dep_dw - dep_up)) 385 ENDIF 386 387 IF ( isodd(je) ) THEN 388 IF ( gst_envl(je) == 0 ) THEN 389 DO jk = 1, num_s 390 ind = s_1st+jk-1 391 ! W-GRID 392 z_gsigw1(ind) = -s_coord(sigma(1,jk,num_s), coefa, & 393 coefb, alpha, num_s, nn_strt(je)) 394 ! T-GRID 395 z_gsigt1(ind) = -s_coord(sigma(0,jk,num_s), coefa, & 396 coefb, alpha, num_s, nn_strt(je)) 397 zcoefw = z_gsigw1(ind) 398 zcoeft = z_gsigt1(ind) 399 400 IF (nn_strt(je) /= 2) THEN 401 gdept_1d(ind) = z_mes(dep_up, dep_dw, zcoeft, & 402 -sigma(0,jk,num_s), env_hc) 403 gdepw_1d(ind) = z_mes(dep_up, dep_dw, zcoefw, & 404 -sigma(1,jk,num_s), env_hc) 405 ELSE 406 gdept_1d(ind) = z_mes(dep_up, dep_dw, zcoeft, & 407 -sigma(0,jk,num_s), 0._wp) 408 gdepw_1d(ind) = z_mes(dep_up, dep_dw, zcoefw, & 409 -sigma(1,jk,num_s), 0._wp) 410 END IF 411 ENDDO 412 ENDIF 413 414 ELSE 415 416 ! Ensure that for splines we do not use SF12 417 !IF (nn_strt(je) > 1) nn_strt(je) = 1 418 419 ! Interval's endpoints TAU are 0 and 1. 420 tau(1) = 0._wp 421 tau(n) = 1._wp 422 423 IF ( je == 2 ) THEN 424 dep0 = 0._wp 425 ELSE 426 dep0 = rn_ebot_max(je-2) 427 END IF 428 dep1 = dep_up 429 dep2 = dep_dw 430 dep3 = rn_ebot_max(je+1) 431 432 ! Computing derivative analytically 433 ! 434 ! 1. Derivative for upper envelope 435 ibcbeg = 1 436 IF (nn_strt(je-1) == 2) THEN 437 alpha = rn_e_al(je-1) 438 num_s_up = nn_slev(je-1) 439 coefa = rn_e_th(je-1) / (dep1-dep0) 440 coefb = rn_e_bb(je-1) + ((dep1-dep0)*rn_e_ba(je-1)) 441 coefb = 1.0_wp-(coefb/(dep1-dep0)) 442 ELSE 443 alpha = 0._wp 444 num_s_up = 1 445 coefa = rn_e_th(je-1) 446 coefb = rn_e_bb(je-1) 447 END IF 448 D1s_up = D1s_coord(-1._wp, coefa, coefb, & 449 alpha, num_s_up, nn_strt(je-1)) 450 d(1,1) = dep_up 451 IF (nn_strt(je-1) == 2) THEN 452 d(2,1) = D1z_mes(dep0, dep1, D1s_up, 0._wp, 1) 453 ELSE 454 d(2,1) = D1z_mes(dep0, dep1, D1s_up, rn_e_hc(je-1), 1) 455 END IF 456 ! 457 ! 2. Derivative for deeper envelope 458 IF ( gst_envl(je+1) == 0 ) THEN 459 ibcend = 1 ! Bound. cond for the 1st derivative 460 IF (nn_strt(je+1) == 2) THEN 461 alpha = rn_e_al(je+1) 462 num_s_dw = nn_slev(je+1) 463 coefa = rn_e_th(je+1) / (dep3-dep2) 464 coefb = rn_e_bb(je+1) + ((dep3-dep2)*rn_e_ba(je+1)) 465 coefb = 1.0_wp-(coefb/(dep3-dep2)) 466 ELSE 467 alpha = 0._wp 468 num_s_dw = 1 469 coefa = rn_e_th(je+1) 470 coefb = rn_e_bb(je+1) 471 END IF 472 D1s_dw = D1s_coord(0._wp, coefa, coefb, & 473 alpha, num_s_dw, nn_strt(je+1)) 474 d(1,2) = dep_dw 475 IF (nn_strt(je+1) == 2) THEN 476 d(2,2) = D1z_mes(dep2, dep3, D1s_dw, 0._wp, 1) 477 ELSE 478 d(2,2) = D1z_mes(dep2, dep3, D1s_dw, rn_e_hc(je+1), 1) 479 END IF 480 ELSE 481 ibcend = 2 ! Bound. cond for the 2nd derivative 482 d(2,2) = 0._wp ! 2nd der = 0 483 ENDIF 484 ! 485 ! Computing spline's coefficients 486 IF ( lwp ) THEN 487 WRITE(numout,*) "BOUNDARY COND. to CONSTRAIN SPLINES:" 488 WRITE(numout,*) "ibcbeg: ", ibcbeg 489 WRITE(numout,*) "ibcend: ", ibcend 490 ENDIF 491 c = cub_spl(tau, d, n, ibcbeg, ibcend ) 492 493 IF ( lwp ) THEN 494 WRITE(numout,*) "Subzone between depth ", dep1, " and depth ", dep2 495 WRITE(numout,*) "" 496 WRITE(numout,*) "D1s_up = ", D1s_up 497 WRITE(numout,*) "1st deriv. of z(mes) at depth ", dep1, " is ", d(2,1) 498 WRITE(numout,*) "D1s_dw = ", D1s_dw 499 WRITE(numout,*) "1st deriv. of z(mes) at depth ", dep2, " is ", d(2,2) 500 WRITE(numout,*) "" 501 WRITE(numout,*) "CUBIC SPLINE COEFFICENTS" 502 WRITE(numout,*) "c1 = ", c(1,1) 503 WRITE(numout,*) "c2 = ", c(2,1) 504 WRITE(numout,*) "c3 = ", c(3,1) 505 WRITE(numout,*) "c4 = ", c(4,1) 506 WRITE(numout,*) "" 507 WRITE(numout,*) "CUB. SPL. at s = ", -sigma(1, 1, num_s), & 508 " is ", z_cub_spl(-sigma(1, 1, num_s), tau, c) 509 WRITE(numout,*) "1st deriv. of CUB. SPL. at s = ", -sigma(1, 1, num_s) , & 510 " is ", D1z_cub_spl(-sigma(1, 1, num_s), tau, c, 1) 511 WRITE(numout,*) "CUB. SPL. at s = ", -sigma(1, num_s, num_s) , & 512 " is ", z_cub_spl(-sigma(1, num_s, num_s), tau, c) 513 WRITE(numout,*) "1st deriv. of CUB. SPL. at s = ", -sigma(1, num_s, num_s) , & 514 " is ", D1z_cub_spl(-sigma(1, num_s, num_s), tau, c, 1) 515 END IF 516 517 env_hc = rn_e_hc(je) 518 coefa = rn_e_th(je) 519 coefb = rn_e_bb(je) 520 alpha = rn_e_al(je) 521 522 DO jk = 1, num_s 523 ind = s_1st+jk-1 524 ! W-GRID 525 z_gsigw1(s_1st+jk-1) = -s_coord( sigma(1,jk,num_s), coefa, & 526 coefb, alpha, num_s, nn_strt(je) ) 527 ! T-GRID 528 z_gsigt1(s_1st+jk-1) = -s_coord( sigma(0,jk,num_s), coefa, & 529 coefb, alpha, num_s, nn_strt(je) ) 530 st = z_gsigt1(ind) 531 sw = z_gsigw1(ind) 532 gdept_1d(ind) = z_cub_spl(st, tau, c) 533 gdepw_1d(ind) = z_cub_spl(sw, tau, c) 534 END DO 535 536 END IF 537 538 END DO 539 540 !================================================================== 541 ! 2. Computing derivative of MEs-coordinate and scale factors 542 ! as finite differences 543 !================================================================== 544 545 ! Computing derivative of MEs-coordinate as finite differences 546 DO jk = 1, jpkm1 547 z_esigt1(jk) = z_gsigw1(jk+1) - z_gsigw1(jk) 548 z_esigw1(jk+1) = z_gsigt1(jk+1) - z_gsigt1(jk) 549 END DO 550 ! Surface 551 z_esigw1( 1 ) = 2._wp * ( z_gsigt1( 1 ) - z_gsigw1( 1 ) ) 552 ! Bottom 553 z_esigt1(jpk) = 2._wp * ( z_gsigt1(jpk) - z_gsigw1(jpk) ) 554 555 ! Computing scale factors as finite differences 556 DO jk=1,jpkm1 557 e3t_1d(jk) = gdepw_1d(jk+1) - gdepw_1d(jk) 558 e3w_1d(jk+1) = gdept_1d(jk+1) - gdept_1d(jk) 559 ENDDO 560 ! Bottom: 561 jk = jpk 562 e3t_1d(jk) = 2.0_wp * (gdept_1d(jk) - gdepw_1d(jk)) 563 ! and Surface (H. Liu) 564 jk = 1 565 e3w_1d(jk) = 2.0_wp * (gdept_1d(1) - gdepw_1d(1)) 566 567 568 ! Control print 569 IF ( lwp ) THEN 570 WRITE(numout,*) "" 571 WRITE(numout,*) "1a. MEs-levels and first derivative" 572 WRITE(numout,*) "" 573 WRITE(numout,*) " k z_gsigw1 z_esigw1 z_gsigt1 z_esigt1" 574 WRITE(numout,*) "" 575 DO jk = 1, jpk 576 WRITE(numout,*) jk, z_gsigw1(jk), z_esigw1(jk), z_gsigt1(jk), z_esigt1(jk) 577 ENDDO 578 WRITE(numout,*) "-----------------------------------------------------------------" 579 580 WRITE(numout,*) "" 581 WRITE(numout,*) "2a. Depth of MEs-levels and scale factors" 582 WRITE(numout,*) "" 583 WRITE(numout,*) " k gdepw_1d e3w_1d gdepw_1d e3w_1d" 584 WRITE(numout,*) "" 585 DO jk = 1, jpk 586 WRITE(numout,*) jk, gdepw_1d(jk), e3w_1d(jk), gdept_1d(jk), e3t_1d(jk) 587 ENDDO 588 WRITE(numout,*) "-----------------------------------------------------------------" 589 ENDIF 590 591 ! Checking monotonicity for cubic splines 592 DO jk = 1, jpk-1 593 IF ( gdept_1d(jk+1) < gdept_1d(jk) ) THEN 594 WRITE(ctlmes,*) 'NOT MONOTONIC gdept_1d: change envelopes' 595 CALL ctl_stop( ctlmes ) 596 END IF 597 IF ( gdepw_1d(jk+1) < gdepw_1d(jk) ) THEN 598 WRITE(ctlmes,*) 'NOT MONOTONIC gdepw_1d: change envelopes' 599 CALL ctl_stop( ctlmes ) 600 END IF 601 END DO 602 603 !======================== 604 ! 3. Compute 3D MEs mesh 605 !======================== 606 607 DO je = 1, tot_env ! LOOP over all the requested envelopes 608 609 IF (je == 1) THEN 610 s_1st = 1 611 num_s = nn_slev(je) 612 env_up = scosrf ! surface 613 ELSE 614 s_1st = s_1st + num_s - 1 615 num_s = nn_slev(je) + 1 616 env_up = envlt(:,:,je-1) 617 ENDIF 618 619 env_dw = envlt(:,:,je) 620 env_hc = rn_e_hc(je) 621 alpha = rn_e_al(je) 622 623 IF ( isodd(je) ) THEN 624 625 IF ( gst_envl(je) == 0 ) THEN 626 DO jk = 1, num_s 627 ind = s_1st+jk-1 628 DO jj = 1,jpj 629 DO ji = 1,jpi 630 IF (env_dw(ji,jj) < env_hc) THEN 631 ! shallow water, uniform sigma 632 zcoefw = -sigma(1, jk, num_s) 633 zcoeft = -sigma(0, jk, num_s) 634 !IF (nn_strt(je) == 2) THEN 635 ! ! shallow water, z coordinates 636 ! ! SF12 in AMM* configurations uses 637 ! ! this (ln_sigcrit=.false.) 638 ! zcoefw = zcoefw*(env_hc/(env_dw(ji,jj)-env_up(ji,jj))) 639 ! zcoeft = zcoeft*(env_hc/(env_dw(ji,jj)-env_up(ji,jj))) 640 !ENDIF 641 ELSE ! deep water, stretched s 642 IF (nn_strt(je) == 2) THEN 643 coefa = rn_e_th(je) / (env_dw(ji,jj)-env_up(ji,jj)) 644 coefb = rn_e_bb(je) + ((env_dw(ji,jj)-env_up(ji,jj))*rn_e_ba(je)) 645 coefb = 1.0_wp-(coefb/(env_dw(ji,jj)-env_up(ji,jj))) 646 ! W-GRID 647 zcoefw = -s_coord( sigma(1,jk,num_s), coefa, & 648 coefb, alpha, num_s, nn_strt(je) ) 649 ! T-GRID 650 zcoeft = -s_coord( sigma(0,jk,num_s), coefa, & 651 coefb, alpha, num_s, nn_strt(je) ) 652 ELSE 653 zcoefw = z_gsigw1(ind) 654 zcoeft = z_gsigt1(ind) 655 END IF 656 ENDIF 657 658 IF (nn_strt(je) /= 2) THEN 659 gdept_0(ji,jj,ind) = z_mes(env_up(ji,jj), env_dw(ji,jj), zcoeft, & 660 -sigma(0, jk, num_s), env_hc) 661 662 gdepw_0(ji,jj,ind) = z_mes(env_up(ji,jj), env_dw(ji,jj), zcoefw, & 663 -sigma(1, jk, num_s), env_hc) 664 ELSE 665 gdept_0(ji,jj,ind) = z_mes(env_up(ji,jj), env_dw(ji,jj), zcoeft, & 666 -sigma(0, jk, num_s), 0._wp) 667 668 gdepw_0(ji,jj,ind) = z_mes(env_up(ji,jj), env_dw(ji,jj), zcoefw, & 669 -sigma(1, jk, num_s), 0._wp) 670 END IF 671 END DO 672 END DO 673 END DO 674 END IF 675 676 ELSE 677 678 ! Interval's endpoints TAU are 0 and 1. 679 tau(1) = 0._wp 680 tau(n) = 1._wp 681 682 IF ( je == 2 ) THEN 683 env0 = scosrf 684 ELSE 685 env0 = envlt(:,:,je-2) 686 END IF 687 env1 = env_up 688 env2 = env_dw 689 env3 = envlt(:,:,je+1) 690 691 DO jj = 1,jpj 692 DO ji = 1,jpi 693 d(1,1) = env_up(ji,jj) 694 d(1,2) = env_dw(ji,jj) 695 ! 696 ! Computing derivative analytically 697 ! 698 ! 1. Derivative for upper envelope 699 ibcbeg = 1 700 IF (nn_strt(je-1) == 2) THEN 701 alpha = rn_e_al(je-1) 702 num_s_up = nn_slev(je-1) 703 coefa = rn_e_th(je-1) / & 704 (env1(ji,jj)-env0(ji,jj)) 705 coefb = rn_e_bb(je-1) + & 706 ((env1(ji,jj)-env0(ji,jj))*rn_e_ba(je-1)) 707 coefb = 1.0_wp-(coefb/(env1(ji,jj)-env0(ji,jj))) 708 ELSE 709 alpha = 0._wp 710 num_s_up = 1 711 coefa = rn_e_th(je-1) 712 coefb = rn_e_bb(je-1) 713 END IF 714 D1s_up = D1s_coord(-1._wp, coefa, coefb, & 715 alpha, num_s_up, nn_strt(je-1)) 716 IF (nn_strt(je-1) == 2) THEN 717 d(2,1) = D1z_mes(env0(ji,jj), env1(ji,jj), D1s_up, & 718 0._wp, 1) 719 ELSE 720 d(2,1) = D1z_mes(env0(ji,jj), env1(ji,jj), D1s_up, & 721 rn_e_hc(je-1), 1) 722 END IF 723 ! 724 ! 2. Derivative for lower envelope 725 IF ( gst_envl(je+1) == 0 ) THEN 726 ibcend = 1 ! Bound. cond for the 1st derivative 727 IF (nn_strt(je+1) == 2) THEN 728 alpha = rn_e_al(je+1) 729 num_s_dw = nn_slev(je+1) 730 coefa = rn_e_th(je+1) / & 731 (env3(ji,jj)-env2(ji,jj)) 732 coefb = rn_e_bb(je+1) + & 733 ((env3(ji,jj)-env2(ji,jj))*rn_e_ba(je+1)) 734 coefb = 1.0_wp-(coefb/(env3(ji,jj)-env2(ji,jj))) 735 ELSE 736 alpha = 0._wp 737 num_s_dw = 1 738 coefa = rn_e_th(je+1) 739 coefb = rn_e_bb(je+1) 740 END IF 741 D1s_dw = D1s_coord(0._wp, coefa, coefb, & 742 alpha, num_s_dw, nn_strt(je+1)) 743 744 IF (nn_strt(je+1) == 2) THEN 745 d(2,2) = D1z_mes(env2(ji,jj), env3(ji,jj), D1s_dw, & 746 0._wp, 1) 747 ELSE 748 d(2,2) = D1z_mes(env2(ji,jj), env3(ji,jj), D1s_dw, & 749 rn_e_hc(je+1), 1) 750 END IF 751 ELSE 752 ibcend = 2 ! Bound. cond for the 2nd derivative 753 d(2,2) = 0._wp ! 2nd der = 0 754 ENDIF 755 756 ! Computing spline's coefficients 757 !IF ( lwp ) THEN 758 ! WRITE(numout,*) "BOUNDARY COND. to CONSTRAIN SPLINES:" 759 ! WRITE(numout,*) "ibcbeg: ", ibcbeg 760 ! WRITE(numout,*) "ibcend: ", ibcend 761 !ENDIF 762 c = cub_spl(tau, d, n, ibcbeg, ibcend ) 763 764 ! Computing model level depths 765 DO jk = 1, num_s 766 ind = s_1st+jk-1 767 st = z_gsigt1(ind) 768 sw = z_gsigw1(ind) 769 gdept_0(ji,jj,ind) = z_cub_spl(st, tau, c) 770 gdepw_0(ji,jj,ind) = z_cub_spl(sw, tau, c) 771 END DO 772 773 ! Control print 774 IF ( lwp ) THEN 775 IF ( jj == 25 .AND. ji == 7 ) THEN 776 WRITE(numout,*) "" 777 ! WRITE(numout,*) "env_up = ", env1(ji,jj) 778 WRITE(numout,*) "1st deriv. at TAU = 0._wp is ", d(2,1) 779 ! WRITE(numout,*) "env_dw = ", env2(ji,jj) 780 WRITE(numout,*) "1st deriv. at TAU = 1._wp is ", d(2,2) 781 ! WRITE(numout,*) "CUBIC SPLINE at s = ", -sigma(1, 1, num_s) , " is ", & 782 ! z_cub_spl(-sigma(1, 1, num_s), tau, c) 783 WRITE(numout,*) "1st deriv. at s = ", -sigma(1, 1, num_s) , " is ", & 784 D1z_cub_spl(-sigma(1, 1, num_s), tau, c, 1) 785 ! WRITE(numout,*) "CUBIC SPLINE at s = ", -sigma(1, num_s, num_s) , " is ", & 786 ! z_cub_spl(-sigma(1, num_s, num_s), tau, c) 787 WRITE(numout,*) "1st deriv. at s = ", -sigma(1, num_s, num_s) , " is ", & 788 D1z_cub_spl(-sigma(1, num_s, num_s), tau, c, 1) 789 END IF 790 END IF 791 792 END DO 793 END DO 794 795 END IF 796 797 END DO 798 799 ! ============================================= 800 ! 4. COMPUTE e3t, e3w for ALL general levels 801 ! as finite differences 802 ! ============================================= 803 ! 804 DO jk=1,jpkm1 805 e3t_0(:,:,jk) = gdepw_0(:,:,jk+1) - gdepw_0(:,:,jk) 806 e3w_0(:,:,jk+1) = gdept_0(:,:,jk+1) - gdept_0(:,:,jk) 807 ENDDO 808 ! Surface 809 jk = 1 810 e3w_0(:,:,jk) = 2.0_wp * (gdept_0(:,:,1) - gdepw_0(:,:,1)) 811 ! 812 ! Bottom 813 jk = jpk 814 e3t_0(:,:,jk) = 2.0_wp * (gdept_0(:,:,jk) - gdepw_0(:,:,jk)) 815 816 IF ( lwp ) THEN 817 jj = 25 818 ji = 7 819 WRITE(numout,*) "" 820 WRITE(numout,*) " k gdepw_0(k) gdept_0(k) e3w_0(k) e3t_0(k)" 821 WRITE(numout,*) "" 822 DO jk = 1, jpk 823 WRITE(numout,*) jk, gdepw_0(ji,jj,jk), gdept_0(ji,jj,jk), e3w_0(ji,jj,jk), e3t_0(ji,jj,jk) 824 ENDDO 825 END IF 826 827 !============================== 828 ! Computing gdept3w 829 !============================== 830 831 gde3w_0(:,:,1) = 0.5 * e3w_0(:,:,1) 832 DO jk = 2, jpk 833 gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk) 834 END DO 835 e3t_0(:,:,jpk) = 2.0_wp * (gdept_0(:,:,jpk) - gdepw_0(:,:,jpk)) 836 837 !============================== 838 ! Computing mbathy 839 !============================== 840 841 IF( lwp ) WRITE(numout,*) '' 842 IF( lwp ) WRITE(numout,*) 'MES mbathy:' 843 IF( lwp ) WRITE(numout,*) '' 844 845 DO jj = 1, jpj 846 DO ji = 1, jpi 847 DO jk = 1, jpkm1 848 !IF( lwp ) WRITE(numout,*) 'scobot: ', scobot(ji,jj), 'gdept_0: ', gdept_0(ji,jj,jk) 849 IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) 850 IF( scobot(ji,jj) == 0.e0 ) mbathy(ji,jj) = 0 851 IF( scobot(ji,jj) < 0.e0 ) mbathy(ji,jj) = int( scobot(ji,jj)) !???? do we need it? 852 END DO 853 !IF( lwp ) WRITE(numout,*) 'mbathy(',ji,',',jj,') = ', mbathy(ji,jj) 854 END DO 855 END DO 856 857 !============================================== 858 ! Computing e3u_0, e3v_0, e3f_0, e3uw_0, e3vw_0 859 !============================================== 860 861 DO jj = 1, jpjm1 862 DO ji = 1, jpim1 863 DO jk = 1, jpk 864 865 e3u_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3t_0(ji,jj,jk) + & 866 MIN(1,mbathy(ji+1,jj))*e3t_0(ji+1,jj,jk) ) / & 867 MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji+1,jj)) ) 868 869 e3v_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3t_0(ji,jj,jk) + & 870 MIN(1,mbathy(ji,jj+1))*e3t_0(ji,jj+1,jk) ) / & 871 MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1)) ) 872 873 e3uw_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3w_0(ji,jj,jk) + & 874 MIN(1,mbathy(ji+1,jj))*e3w_0(ji+1,jj,jk) ) / & 875 MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji+1,jj)) ) 876 877 e3vw_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3w_0(ji,jj,jk) + & 878 MIN(1,mbathy(ji,jj+1))*e3w_0(ji,jj+1,jk) ) / & 879 MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1)) ) 880 881 e3f_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3t_0(ji,jj,jk) + & 882 MIN(1,mbathy(ji+1,jj))*e3t_0(ji+1,jj,jk) + & 883 MIN(1,mbathy(ji+1,jj+1))* e3t_0(ji+1,jj+1,jk) + & 884 MIN(1,mbathy(ji,jj+1))*e3t_0(ji,jj+1,jk) ) / & 885 MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1)) & 886 + MIN(1,mbathy(ji+1,jj))+MIN(1,mbathy(ji+1,jj+1)) ) 887 ENDDO 888 ENDDO 889 ENDDO 890 891 ! Adjusting for geopotential levels, if any 892 893 DO je = 1, tot_env ! LOOP over all the requested envelopes 894 895 IF (je == 1) THEN 896 s_1st = 1 897 num_s = nn_slev(je) 898 max_env_up = 0._wp ! surface 899 min_env_up = 0._wp ! surface 900 ELSE 901 s_1st = s_1st + num_s - 1 902 num_s = nn_slev(je) + 1 903 max_env_up = rn_ebot_max(je-1) 904 min_env_up = rn_ebot_min(je-1) 905 ENDIF 906 907 max_env_dw = rn_ebot_max(je) 908 min_env_dw = rn_ebot_min(je) 909 910 IF ( max_env_up == min_env_up .AND. max_env_dw == min_env_dw ) THEN 911 912 DO jj = 1,jpjm1 913 DO ji = 1,jpim1 914 DO jk = 1, num_s 915 916 ind = s_1st+jk-1 917 e3u_0 (ji,jj,ind) = MIN( e3t_0(ji,jj,ind), e3t_0(ji+1,jj,ind)) 918 e3uw_0(ji,jj,ind) = MIN( e3w_0(ji,jj,ind), e3w_0(ji+1,jj,ind)) 919 e3v_0 (ji,jj,ind) = MIN( e3t_0(ji,jj,ind), e3t_0(ji,jj+1,ind)) 920 e3vw_0(ji,jj,ind) = MIN( e3w_0(ji,jj,ind), e3w_0(ji,jj+1,ind)) 921 e3f_0 (ji,jj,ind) = MIN( e3t_0(ji,jj,ind), e3t_0(ji+1,jj,ind)) 922 923 ENDDO 924 ENDDO 925 ENDDO 926 927 ENDIF 928 929 ENDDO 930 931 !============================== 932 ! Computing gdept3w 933 !============================== 934 935 gde3w_0(:,:,1) = 0.5 * e3w_0(:,:,1) 936 DO jk = 2, jpk 937 gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk) 938 END DO 939 940 !================================================================================= 941 ! From here equal to sco code - domzgr.F90 line 2079 942 ! Lateral B.C. 943 944 CALL lbc_lnk( e3t_0 , 'T', 1._wp ) 945 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) 946 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) 947 CALL lbc_lnk( e3f_0 , 'F', 1._wp ) 948 CALL lbc_lnk( e3w_0 , 'W', 1._wp ) 949 CALL lbc_lnk( e3uw_0, 'U', 1._wp ) 950 CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 951 952 where (e3t_0 (:,:,:).eq.0.0) e3t_0(:,:,:) = 1.0 953 where (e3u_0 (:,:,:).eq.0.0) e3u_0(:,:,:) = 1.0 954 where (e3v_0 (:,:,:).eq.0.0) e3v_0(:,:,:) = 1.0 955 where (e3f_0 (:,:,:).eq.0.0) e3f_0(:,:,:) = 1.0 956 where (e3w_0 (:,:,:).eq.0.0) e3w_0(:,:,:) = 1.0 957 where (e3uw_0 (:,:,:).eq.0.0) e3uw_0(:,:,:) = 1.0 958 where (e3vw_0 (:,:,:).eq.0.0) e3vw_0(:,:,:) = 1.0 959 960 IF( lwp ) WRITE(numout,*) 'Refine mbathy' 961 DO jj = 1, jpj 962 DO ji = 1, jpi 963 DO jk = 1, jpkm1 964 IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) 965 END DO 966 IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0 967 END DO 968 END DO 969 970 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) ' MIN val mbathy h90 ', MINVAL( mbathy(:,:) ), & 971 & ' MAX ', MAXVAL( mbathy(:,:) ) 972 973 IF( nprint == 1 .AND. lwp ) THEN ! min max values over the local domain 974 WRITE(numout,*) 'zgr_sco : mbathy' 975 WRITE(numout,*) '~~~~~~~' 976 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 977 WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), & 978 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ' , MINVAL( gde3w_0(:,:,:) ) 979 WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 (:,:,:) ), & 980 & ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 (:,:,:) ), & 981 & ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 (:,:,:) ), & 982 & ' w ', MINVAL( e3w_0 (:,:,:) ) 983 984 WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), & 985 & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ' , MAXVAL( gde3w_0(:,:,:) ) 986 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 (:,:,:) ), & 987 & ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 (:,:,:) ), & 988 & ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 (:,:,:) ), & 989 & ' w ', MAXVAL( e3w_0 (:,:,:) ) 990 ENDIF 991 992 !================================================================================ 993 ! Check coordinates makes sense 994 !================================================================================ 995 996 DO ji = 1, jpi 997 DO jj = 1, jpj 998 DO jk = 1, mbathy(ji,jj) 999 ! check coordinate is monotonically increasing 1000 IF (e3w_0(ji,jj,jk) <= 0._wp .OR. e3t_0(ji,jj,jk) <= 0._wp ) THEN 1001 WRITE(ctmp1,*) 'ERROR zgr_mes : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 1002 WRITE(numout,*) 'ERROR zgr_mes : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 1003 WRITE(numout,*) 'e3w',e3w_0(ji,jj,:) 1004 WRITE(numout,*) 'e3t',e3t_0(ji,jj,:) 1005 CALL ctl_stop( ctmp1 ) 1006 ENDIF 1007 ! and check it has never gone negative 1008 IF ( gdepw_0(ji,jj,jk) < 0._wp .OR. gdept_0(ji,jj,jk) < 0._wp ) THEN 1009 WRITE(ctmp1,*) 'ERROR zgr_mes : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 1010 WRITE(numout,*) 'ERROR zgr_mes : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 1011 WRITE(numout,*) 'gdepw',gdepw_0(ji,jj,:) 1012 WRITE(numout,*) 'gdept',gdept_0(ji,jj,:) 1013 CALL ctl_stop( ctmp1 ) 1014 ENDIF 1015 ! ! and check it never exceeds the total depth 1016 ! IF ( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN 1017 ! WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk 1018 ! WRITE(numout,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk 1019 ! WRITE(numout,*) 'gdepw',fsdepw(ji,jj,:) 1020 ! CALL ctl_stop( ctmp1 ) 1021 ! ENDIF 1022 END DO 1023 1024 ! DO jk = 1, mbathy(ji,jj)-1 1025 ! ! and check it never exceeds the total depth 1026 ! IF ( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN 1027 ! WRITE(ctmp1,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk 1028 ! WRITE(numout,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk 1029 ! WRITE(numout,*) 'gdept',fsdept(ji,jj,:) 1030 ! CALL ctl_stop( ctmp1 ) 1031 ! ENDIF 1032 ! END DO 1033 END DO 1034 END DO 1035 ! 1036 CALL wrk_dealloc( jpi, jpj, max_nn_env , envlt ) 1037 CALL wrk_dealloc( jpi, jpj, env_up, env_dw ) 1038 CALL wrk_dealloc( jpk, z_gsigw1, z_gsigt1, z_esigw1, z_esigt1 ) 1039 ! 1040 IF( nn_timing == 1 ) CALL timing_stop('zgr_mes') 1041 1042 END SUBROUTINE zgr_mes_build 1043 1044 ! ===================================================================================================== 1045 1046 SUBROUTINE zgr_dom_read(ld_isfcav, & ! ice-cavities 1047 & pdept_1d, pdepw_1d, pe3t_1d , pe3w_1d, & ! 1D reference vertical coordinate 1048 & pdept , pdepw , & ! 3D t & w-points depth 1049 & pe3t , pe3u , pe3v , pe3f , & ! vertical scale factors 1050 & pe3w , pe3uw , pe3vw ) ! - - - 84 SUBROUTINE zgr_dom_read 1051 85 !!--------------------------------------------------------------------- 1052 86 !! *** ROUTINE zgr_dom_read *** … … 1055 89 !! 1056 90 !!---------------------------------------------------------------------- 1057 LOGICAL , INTENT(out) :: ld_isfcav ! under iceshelf cavity flag1058 REAL(wp), DIMENSION(:) , INTENT(out) :: pdept_1d, pdepw_1d ! 1D grid-point depth [m]1059 REAL(wp), DIMENSION(:) , INTENT(out) :: pe3t_1d , pe3w_1d ! 1D vertical scale factors [m]1060 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pdept, pdepw ! grid-point depth [m]1061 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pe3t , pe3u , pe3v , pe3f ! vertical scale factors [m]1062 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pe3w , pe3uw, pe3vw ! - - -1063 !1064 91 INTEGER :: jk ! dummy loop index 1065 92 INTEGER :: inum ! local logical unit … … 1078 105 ! !* ocean cavities under iceshelves 1079 106 CALL iom_get( inum, 'ln_isfcav', z_cav ) 1080 IF( z_cav == 0._wp ) THEN ; l d_isfcav = .false. ; ELSE ; ld_isfcav = .true. ; ENDIF107 IF( z_cav == 0._wp ) THEN ; ln_isfcav = .false. ; ELSE ; ln_isfcav = .true. ; ENDIF 1081 108 ! 1082 109 ! !* vertical scale factors 1083 CALL iom_get( inum, jpdom_unknown, 'e3t_1d' , pe3t_1d ) ! 1D reference coordinate 1084 CALL iom_get( inum, jpdom_unknown, 'e3w_1d' , pe3w_1d ) 1085 ! 1086 CALL iom_get( inum, jpdom_global, 'e3t_0' , pe3t ) ! 3D coordinate 1087 CALL iom_get( inum, jpdom_global, 'e3u_0' , pe3u ) 1088 CALL iom_get( inum, jpdom_global, 'e3v_0' , pe3v ) 1089 CALL iom_get( inum, jpdom_global, 'e3f_0' , pe3f ) 1090 CALL iom_get( inum, jpdom_global, 'e3w_0' , pe3w ) 1091 CALL iom_get( inum, jpdom_global, 'e3uw_0' , pe3uw ) 1092 CALL iom_get( inum, jpdom_global, 'e3vw_0' , pe3vw ) 110 CALL iom_get( inum, jpdom_unknown, 'e3t_1d' , e3t_1d_in ) ! 1D reference coordinate 111 CALL iom_get( inum, jpdom_unknown, 'e3w_1d' , e3w_1d_in ) 112 ! 113 CALL iom_get( inum, jpdom_global, 'bottom_level' , mbathy_in ) ! 2D mbathy 114 ! 115 CALL iom_get( inum, jpdom_global, 'e3t_0' , e3t_in ) ! 3D coordinate 116 CALL iom_get( inum, jpdom_global, 'e3u_0' , e3u_in ) 117 CALL iom_get( inum, jpdom_global, 'e3v_0' , e3v_in ) 118 CALL iom_get( inum, jpdom_global, 'e3f_0' , e3f_in ) 119 CALL iom_get( inum, jpdom_global, 'e3w_0' , e3w_in ) 120 CALL iom_get( inum, jpdom_global, 'e3uw_0' , e3uw_in ) 121 CALL iom_get( inum, jpdom_global, 'e3vw_0' , e3vw_in ) 1093 122 ! 1094 123 ! !* depths … … 1100 129 CALL ctl_warn( 'zgr_dom_read : old definition of depths and scale factors used ', & 1101 130 & ' depths at t- and w-points read in the domain configuration file') 1102 CALL iom_get( inum, jpdom_unknown, 'gdept_1d', pdept_1d)1103 CALL iom_get( inum, jpdom_unknown, 'gdepw_1d', pdepw_1d)1104 CALL iom_get( inum, jpdom_global , 'gdept_0' , pdept)1105 CALL iom_get( inum, jpdom_global , 'gdepw_0' , pdepw)131 CALL iom_get( inum, jpdom_unknown, 'gdept_1d', gdept_1d_in ) 132 CALL iom_get( inum, jpdom_unknown, 'gdepw_1d', gdepw_1d_in ) 133 CALL iom_get( inum, jpdom_global , 'gdept_0' , gdept_in ) 134 CALL iom_get( inum, jpdom_global , 'gdepw_0' , gdepw_in ) 1106 135 ! 1107 136 ELSE !- depths computed from e3. scale factors 1108 CALL e3_to_depth( pe3t_1d, pe3w_1d, pdept_1d, pdepw_1d) ! 1D reference depth1109 CALL e3_to_depth( pe3t , pe3w , pdept , pdepw) ! 3D depths137 CALL e3_to_depth( e3t_1d_in, e3w_1d_in, gdept_1d_in, gdepw_1d_in ) ! 1D reference depth 138 CALL e3_to_depth( e3t_in , e3w_in , gdept_in , gdepw_in ) ! 3D depths 1110 139 IF(lwp) THEN 1111 140 WRITE(numout,*) 1112 141 WRITE(numout,*) ' GLOBAL reference 1D z-coordinate depth and scale factors:' 1113 142 WRITE(numout, "(9x,' level gdept_1d gdepw_1d e3t_1d e3w_1d ')" ) 1114 WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, pdept_1d(jk), pdepw_1d(jk), pe3t_1d(jk), pe3w_1d(jk), jk = 1, jpk )143 WRITE(numout, "(10x, i4, 4f9.2)" ) ( jk, gdept_1d_in(jk), gdepw_1d_in(jk), e3t_1d_in(jk), e3w_1d_in(jk), jk = 1, jpk ) 1115 144 ENDIF 1116 145 ENDIF … … 1122 151 ! ===================================================================================================== 1123 152 1124 FUNCTION isodd( a ) RESULT(odd) 1125 ! Notice that this method uses a btest intrinsic function. If you look at any number in binary 1126 ! notation, you'll note that the Least Significant Bit (LSB) will always be set for odd numbers 1127 ! and cleared for even numbers. Hence, all that is necessary to determine if a number is even or 1128 ! odd is to check the LSB. Since bitwise operations like AND, OR, XOR etc. are usually much faster 1129 ! than the arithmetic operations, this method will run a lot faster than the one with modulo 1130 ! operation. 1131 ! http://forums.devshed.com/software-design-43/quick-algorithm-determine-odd-29843.html 1132 1133 INTEGER, INTENT (in) :: a 1134 LOGICAL :: odd 1135 odd = btest(a, 0) 1136 END FUNCTION isodd 1137 1138 ! ===================================================================================================== 1139 1140 FUNCTION iseven( a ) RESULT(even) 1141 INTEGER, INTENT (in) :: a 1142 LOGICAL :: even 1143 even = .NOT. isodd(a) 1144 END FUNCTION iseven 1145 1146 ! ===================================================================================================== 1147 1148 FUNCTION sech( x ) RESULT( seh ) 1149 1150 REAL(wp), INTENT(in ) :: x 1151 REAL(wp) :: seh 1152 1153 seh = 1._wp / COSH(x) 1154 1155 END FUNCTION sech 1156 1157 ! ===================================================================================================== 1158 1159 FUNCTION sigma( vgrid, pk, kmax ) RESULT( ps1 ) 1160 !!---------------------------------------------------------------------- 1161 !! *** ROUTINE sigma *** 1162 !! 1163 !! ** Purpose : provide the analytical function for sigma-coordinate 1164 !! (not stretched s-coordinate). 1165 !! 1166 !! ** Method : the function provide the non-dimensional position of 1167 !! T and W points (i.e. between 0 and 1). 1168 !! T-points at integer values (between 1 and jpk) 1169 !! W-points at integer values - 1/2 (between 0.5 and 1170 !! jpk-0.5) 1171 !! 1172 !!---------------------------------------------------------------------- 1173 INTEGER, INTENT (in) :: pk ! continuous k coordinate 1174 INTEGER, INTENT (in) :: kmax ! number of levels 1175 INTEGER, INTENT (in) :: vgrid ! type of vertical grid: 0 = T, 1 = W 1176 REAL(wp) :: kindx ! index of T or W points 1177 REAL(wp) :: ps1 ! value of sigma coordinate (-1 <= ps1 <=0) 1178 ! 1179 SELECT CASE (vgrid) 1180 1181 CASE (0) ! T-points at integer values (between 1 and jpk) 1182 kindx = REAL(pk,wp) 1183 CASE (1) ! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 1184 kindx = REAL(pk,wp) - 0.5_wp 1185 CASE DEFAULT 1186 WRITE(ctlmes,*) 'No valid vertical grid option selected' 1187 1188 END SELECT 1189 1190 1191 ps1 = -(kindx - 0.5) / REAL(kmax-1) ! The surface is at the first W level 1192 ! while the bottom coincides with the 1193 ! last W level 1194 END FUNCTION sigma 1195 1196 ! ===================================================================================================== 1197 1198 FUNCTION s_coord( s, ca, cb, alpha, kmax, ftype ) RESULT ( pf1 ) 1199 !!---------------------------------------------------------------------- 1200 !! *** ROUTINE s_coord *** 1201 !! 1202 !! ** Purpose : provide the analytical stretching function 1203 !! for s-coordinate. 1204 !! 1205 !! ** Method : if ftype = 2: Siddorn and Furner 2012 1206 !! analytical function is used 1207 !! if ftype = 1: Song and Haidvogel 1994 1208 !! analytical function is used 1209 !! if ftype = 0: Madec et al. 1996 1210 !! analytical function is used 1211 !! (pag 65 of NEMO Manual) 1212 !! Reference : Madec, Lott, Delecluse and 1213 !! Crepon, 1996. JPO, 26, 1393-1408 1214 !! 1215 !! s MUST be NEGATIVE: -1 <= s <= 0 1216 !!---------------------------------------------------------------------- 1217 REAL(wp), INTENT(in) :: s ! continuous not stretched 1218 ! "s" coordinate 1219 REAL(wp), INTENT(in) :: ca ! surface stretch. coeff 1220 REAL(wp), INTENT(in) :: cb ! bottom stretch. coeff 1221 REAL(wp), INTENT(in) :: alpha ! alpha stretch. coeff for SF12 1222 INTEGER, INTENT (in) :: kmax ! number of levels 1223 INTEGER, INTENT(in) :: ftype ! type of stretching function 1224 REAL(wp) :: pf1 ! "s" stretched value 1225 ! SF12 stretch. funct. parameters 1226 REAL(wp) :: psmth ! smoothing parameter for transition 1227 ! between shallow and deep areas 1228 REAL(wp) :: za1,za2,za3 ! local variables 1229 REAL(wp) :: zn1,zn2,ps ! local variables 1230 REAL(wp) :: za,zb,zx ! local variables 1231 !!---------------------------------------------------------------------- 1232 SELECT CASE (ftype) 1233 1234 CASE (0) ! M96 stretching function 1235 pf1 = ( TANH( ca * ( s + cb ) ) - TANH( cb * ca ) ) & 1236 & * ( COSH( ca ) & 1237 & + COSH( ca * ( 2.e0 * cb - 1.e0 ) ) ) & 1238 & / ( 2._wp * SINH( ca ) ) 1239 CASE (1) ! SH94 stretching function 1240 IF ( ca == 0 ) THEN ! uniform sigma 1241 pf1 = s 1242 ELSE ! stretched sigma 1243 pf1 = (1._wp - cb) * SINH(ca*s) / SINH(ca) + & 1244 & cb * ( ( TANH(ca*(s + 0.5_wp)) - TANH(0.5_wp*ca) ) / & 1245 & (2._wp*TANH(0.5_wp*ca)) ) 153 SUBROUTINE zgr_mes_local 154 !!--------------------------------------------------------------------- 155 !! *** ROUTINE zgr_dom_merge *** 156 !! 157 !! ** Purpose : Create a local MEs grid within a gloabal grid 158 !! using different vertical coordinates. 159 !! 160 !!---------------------------------------------------------------------- 161 INTEGER :: ji, jj, jk ! dummy loop index 162 INTEGER :: inum ! local logical unit 163 REAL(wp), DIMENSION(jpi,jpj) :: s2z_msk ! mask for MEs-area (=2), 164 ! transition area (=1), 165 ! z-area (=0) 166 REAL(wp), DIMENSION(jpi,jpj) :: s2z_wgt ! weigths for computing model 167 ! levels in transition area 168 !!---------------------------------------------------------------------- 169 170 IF(lwp) THEN 171 WRITE(numout,*) 172 WRITE(numout,*) ' zgr_mes_merge : read the vertical coordinates in domain_cfg_global.nc file' 173 WRITE(numout,*) ' ~~~~~~~~' 174 ENDIF 175 ! 176 CALL iom_open( 'bathy_meter.nc', inum ) 177 CALL iom_get( inum, jpdom_data, 's2z_msk', s2z_msk) 178 CALL iom_get( inum, jpdom_data, 's2z_wgt', s2z_wgt) 179 180 DO jj = 1,jpj 181 DO ji = 1,jpi 182 SELECT CASE (INT(s2z_msk(ji,jj))) 183 CASE (0) ! global zps area 184 mbathy (ji,jj ) = mbathy_in(ji,jj) 185 gdept_0(ji,jj,:) = gdept_in(ji,jj,:) 186 gdepw_0(ji,jj,:) = gdepw_in(ji,jj,:) 187 e3t_0 (ji,jj,:) = e3t_in (ji,jj,:) 188 e3w_0 (ji,jj,:) = e3w_in (ji,jj,:) 189 e3u_0 (ji,jj,:) = e3u_in (ji,jj,:) 190 e3v_0 (ji,jj,:) = e3v_in (ji,jj,:) 191 e3f_0 (ji,jj,:) = e3f_in (ji,jj,:) 192 e3uw_0 (ji,jj,:) = e3uw_in (ji,jj,:) 193 e3vw_0 (ji,jj,:) = e3vw_in (ji,jj,:) 194 CASE (1) ! MEs to zps transition area 195 gdept_0(ji,jj,:) = s2z_wgt(ji,jj) * gdept_0(ji,jj,:) + & 196 & ( 1._wp - s2z_wgt(ji,jj) ) * gdept_in(ji,jj,:) 197 gdepw_0(ji,jj,:) = s2z_wgt(ji,jj) * gdepw_0(ji,jj,:) + & 198 & ( 1._wp - s2z_wgt(ji,jj) ) * gdepw_in(ji,jj,:) 199 CASE (2) ! MEs area 200 CYCLE 201 END SELECT 202 END DO 203 END DO 204 ! 205 ! e3t, e3w for transition zone 206 ! as finite differences 207 DO jj = 1, jpj 208 DO ji = 1, jpi 209 IF ( s2z_msk(ji,jj) == 1._wp ) THEN 210 DO jk = 1,jpkm1 211 e3t_0(ji,jj,jk) = gdepw_0(ji,jj,jk+1) - gdepw_0(ji,jj,jk) 212 e3w_0(ji,jj,jk+1) = gdept_0(ji,jj,jk+1) - gdept_0(ji,jj,jk) 213 ENDDO 214 ! Surface 215 jk = 1 216 e3w_0(ji,jj,jk) = 2.0_wp * (gdept_0(ji,jj,1) - gdepw_0(ji,jj,1)) 217 ! 218 ! Bottom 219 jk = jpk 220 e3t_0(ji,jj,jk) = 2.0_wp * (gdept_0(ji,jj,jk) - gdepw_0(ji,jj,jk)) 1246 221 END IF 1247 CASE (2) ! SF12 stretching function 1248 psmth = 1.0_wp ! We consider only the case for efold = 0 1249 ps = -s 1250 zn1 = 1. / REAL(kmax-1) 1251 zn2 = 1. - zn1 1252 1253 za1 = (alpha+2.0_wp)*zn1**(alpha+1.0_wp)-(alpha+1.0_wp)*zn1**(alpha+2.0_wp) 1254 za2 = (alpha+2.0_wp)*zn2**(alpha+1.0_wp)-(alpha+1.0_wp)*zn2**(alpha+2.0_wp) 1255 za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1) 1256 1257 za = cb - za3*(ca-za1)-za2 1258 za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp)) ) 1259 zb = (ca - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1) 1260 zx = 1.0_wp-za/2.0_wp-zb 1261 1262 pf1 = za*(ps*(1.0_wp-ps/2.0_wp))+zb*ps**3.0_wp + & 1263 & zx*( (alpha+2.0_wp)*ps**(alpha+1.0_wp) - & 1264 & (alpha+1.0_wp)*ps**(alpha+2.0_wp) ) 1265 pf1 = -pf1*psmth+ps*(1.0_wp-psmth) 1266 CASE (3) ! stretching function for deeper part of the domain 1267 IF ( ca == 0 ) THEN ! uniform sigma 1268 pf1 = s 1269 ELSE ! stretched sigma 1270 pf1 = (1._wp - cb) * SINH(ca*s) / SINH(ca) 222 END DO 223 END DO 224 ! 225 ! MBATHY transition zone 226 DO jj = 1, jpj 227 DO ji = 1, jpi 228 IF ( s2z_msk(ji,jj) == 1._wp ) THEN 229 DO jk = 1, jpkm1 230 IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) 231 IF( scobot(ji,jj) == 0.e0 ) mbathy(ji,jj) = 0 232 IF( scobot(ji,jj) < 0.e0 ) mbathy(ji,jj) = INT( scobot(ji,jj)) ! do we need it? 233 END DO 1271 234 END IF 1272 END SELECT 1273 1274 END FUNCTION s_coord 1275 1276 ! ===================================================================================================== 1277 1278 FUNCTION D1s_coord( s, ca, cb, alpha, kmax, ftype) RESULT( pf1 ) 1279 !!---------------------------------------------------------------------- 1280 !! *** ROUTINE D1s_coord *** 1281 !! 1282 !! ** Purpose : provide the 1st derivative of the analytical 1283 !! stretching function for s-coordinate. 1284 !! 1285 !! ** Method : if ftype = 2: Siddorn and Furner 2012 1286 !! analytical function is used 1287 !! if ftype = 1: Song and Haidvogel 1994 1288 !! analytical function is used 1289 !! if ftype = 0: Madec et al. 1996 1290 !! analytical function is used 1291 !! (pag 65 of NEMO Manual) 1292 !! Reference : Madec, Lott, Delecluse and 1293 !! Crepon, 1996. JPO, 26, 1393-1408 1294 !! 1295 !! s MUST be NEGATIVE: -1 <= s <= 0 1296 !!---------------------------------------------------------------------- 1297 REAL(wp), INTENT(in ) :: s ! continuous not stretched 1298 ! "s" coordinate 1299 REAL(wp), INTENT(in) :: ca ! surface stretch. coeff 1300 REAL(wp), INTENT(in) :: cb ! bottom stretch. coeff 1301 REAL(wp), INTENT(in) :: alpha ! alpha stretch. coeff for SF12 1302 INTEGER, INTENT (in) :: kmax ! number of levels 1303 INTEGER, INTENT(in ) :: ftype ! type of stretching function 1304 REAL(wp) :: pf1 ! "s" stretched value 1305 ! SF12 stretch. funct. parameters 1306 REAL(wp) :: psmth ! smoothing parameter for transition 1307 ! between shallow and deep areas 1308 REAL(wp) :: za1,za2,za3 ! local variables 1309 REAL(wp) :: zn1,zn2,ps ! local variables 1310 REAL(wp) :: za,zb,zx ! local variables 1311 REAL(wp) :: zt1,zt2,zt3 ! local variables 1312 !!---------------------------------------------------------------------- 1313 SELECT CASE (ftype) 1314 1315 CASE (0) ! M96 stretching function 1316 pf1 = ( ca * ( COSH(ca*(2._wp * cb - 1._wp)) + COSH(ca)) * & 1317 sech(ca * (s+cb))**2 ) / (2._wp * SINH(ca)) 1318 CASE (1) ! SH94 stretching function 1319 IF ( ca == 0 ) then ! uniform sigma 1320 pf1 = 1._wp / REAL(kmax,wp) 1321 ELSE ! stretched sigma 1322 pf1 = (1._wp - cb) * ca * COSH(ca*s) / (SINH(ca) * REAL(kmax,wp)) + & 1323 cb * ca * & 1324 (sech((s+0.5_wp)*ca)**2) / (2._wp * TANH(0.5_wp*ca) * REAL(kmax,wp)) 235 END DO 236 END DO 237 ! 238 ! Computing e3u_0, e3v_0, e3f_0, e3uw_0, e3vw_0 239 ! for transition zone 240 ! 241 DO jj = 1, jpjm1 242 DO ji = 1, jpim1 243 IF ( s2z_msk(ji,jj) == 1._wp ) THEN 244 DO jk = 1, jpk 245 e3u_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3t_0(ji,jj,jk) + & 246 MIN(1,mbathy(ji+1,jj))*e3t_0(ji+1,jj,jk) ) / & 247 MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji+1,jj)) ) 248 249 e3v_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3t_0(ji,jj,jk) + & 250 MIN(1,mbathy(ji,jj+1))*e3t_0(ji,jj+1,jk) ) / & 251 MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1)) ) 252 253 e3uw_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3w_0(ji,jj,jk) + & 254 MIN(1,mbathy(ji+1,jj))*e3w_0(ji+1,jj,jk) ) / & 255 MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji+1,jj)) ) 256 257 e3vw_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3w_0(ji,jj,jk) + & 258 MIN(1,mbathy(ji,jj+1))*e3w_0(ji,jj+1,jk) ) / & 259 MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1)) ) 260 261 e3f_0(ji,jj,jk)=(MIN(1,mbathy(ji,jj))* e3t_0(ji,jj,jk) + & 262 MIN(1,mbathy(ji+1,jj))*e3t_0(ji+1,jj,jk) + & 263 MIN(1,mbathy(ji+1,jj+1))* e3t_0(ji+1,jj+1,jk) + & 264 MIN(1,mbathy(ji,jj+1))*e3t_0(ji,jj+1,jk) ) / & 265 MAX( 1, MIN(1,mbathy(ji,jj))+MIN(1,mbathy(ji,jj+1)) & 266 + MIN(1,mbathy(ji+1,jj))+MIN(1,mbathy(ji+1,jj+1)) ) 267 END DO 1325 268 END IF 1326 CASE (2) ! SF12 stretching function 1327 ps = -s 1328 zn1 = 1. / REAL(kmax-1) 1329 zn2 = 1. - zn1 1330 1331 za1 = (alpha+2.0_wp)*zn1**(alpha+1.0_wp)-(alpha+1.0_wp)*zn1**(alpha+2.0_wp) 1332 za2 = (alpha+2.0_wp)*zn2**(alpha+1.0_wp)-(alpha+1.0_wp)*zn2**(alpha+2.0_wp) 1333 za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1) 1334 1335 za = cb - za3*(ca-za1)-za2 1336 za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp)) ) 1337 zb = (ca - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1) 1338 zx = (alpha+2.0_wp)*(alpha+1.0_wp) 1339 1340 zt1 = 0.5_wp*za*(ps-1._wp)*((ps**2.0_wp + 3._wp*ps + 2._wp)*ps**alpha - 2._wp) 1341 zt2 = zb*(zx*ps**(alpha+1.0_wp) - zx*ps**(alpha) + 3._wp*ps**2._wp) 1342 zt3 = zx*(ps-1._wp)*ps**(alpha) 1343 1344 pf1 = zt1 + zt2 - zt3 1345 1346 END SELECT 1347 1348 END FUNCTION D1s_coord 1349 1350 ! ===================================================================================================== 1351 1352 FUNCTION z_mes(dep_up, dep_dw, Cs, s, hc) RESULT( z ) 1353 !!---------------------------------------------------------------------- 1354 !! *** ROUTINE z_mes *** 1355 !! 1356 !! ** Purpose : provide the analytical trasformation from the 1357 !! computational space (mes-coordinate) to the 1358 !! physical space (depth z) 1359 !! 1360 !! N.B.: z is downward positive defined as well as envelope surfaces. 1361 !! Therefore, Cs MUST be positive, meaning 0. <= Cs <= 1. 1362 !!---------------------------------------------------------------------- 1363 REAL(wp), INTENT(in ) :: dep_up ! shallower envelope 1364 REAL(wp), INTENT(in ) :: dep_dw ! deeper envelope 1365 REAL(wp), INTENT(in ) :: Cs ! stretched s coordinate 1366 REAL(wp), INTENT(in ) :: s ! not stretched s coordinate 1367 REAL(wp), INTENT(in ) :: hc ! critical depth 1368 REAL :: z ! downward positive depth 1369 !!---------------------------------------------------------------------- 1370 ! 1371 1372 z = dep_up + hc * s + Cs * (dep_dw - hc - dep_up) 1373 1374 END FUNCTION z_mes 1375 1376 ! ===================================================================================================== 1377 1378 FUNCTION D1z_mes(dep_up, dep_dw, d1Cs, hc, kmax) RESULT( d1z ) 1379 !!---------------------------------------------------------------------- 1380 !! *** ROUTINE D1z_mes *** 1381 !! 1382 !! ** Purpose : provide the 1st derivative of the analytical 1383 !! trasformation from the computational space 1384 !! (mes-coordinate) to the physical space (depth z) 1385 !! 1386 !! N.B.: z is downward positive defined as well as envelope surfaces. 1387 !! Therefore, Cs MUST be positive, meaning 0. <= Cs <= 1. 1388 !!---------------------------------------------------------------------- 1389 REAL(wp), INTENT(in ) :: dep_up ! shallower envelope 1390 REAL(wp), INTENT(in ) :: dep_dw ! deeper envelope 1391 REAL(wp), INTENT(in ) :: d1Cs ! 1st derivative of the 1392 ! stretched s coordinate 1393 REAL(wp), INTENT(in ) :: hc ! critical depth 1394 INTEGER, INTENT(in ) :: kmax ! number of levels 1395 REAL(wp) :: d1z ! 1st derivative of downward 1396 ! positive depth 1397 !!---------------------------------------------------------------------- 1398 ! 1399 1400 d1z = hc / REAL(kmax,wp) + d1Cs * (dep_dw - hc - dep_up) 1401 !d1z = hc + d1Cs * (dep_dw - hc - dep_up) 1402 1403 END FUNCTION D1z_mes 1404 1405 ! ===================================================================================================== 1406 1407 FUNCTION cub_spl(tau, d, n, ibcbeg, ibcend) RESULT ( c ) 1408 !!---------------------------------------------------------------------- 1409 !! 1410 !! CUBSPL defines an interpolatory cubic spline. 1411 !! 1412 !! Discussion: 1413 !! 1414 !! A tridiagonal linear system for the unknown slopes S(I) of 1415 !! F at TAU(I), I=1,..., N, is generated and then solved by Gauss 1416 !! elimination, with S(I) ending up in C(2,I), for all I. 1417 !! 1418 !! Author: Carl de Boor 1419 !! 1420 !! Reference: Carl de Boor, Practical Guide to Splines, 1421 !! Springer, 2001, ISBN: 0387953663, LC: QA1.A647.v27. 1422 !! 1423 !! Parameters: 1424 !! 1425 !! Input: 1426 !! TAU(N) : the abscissas or X values of the data points. 1427 !! The entries of TAU are assumed to be strictly 1428 !! increasing. 1429 !! 1430 !! N : the number of data points. N is assumed to be 1431 !! at least 2. 1432 !! 1433 !! IBCBEG : boundary condition indicator at TAU(1). 1434 !! = 0 no boundary condition at TAU(1) is given. 1435 !! In this case, the "not-a-knot condition" 1436 !! is used. 1437 !! = 1 the 1st derivative at TAU(1) is equal to 1438 !! the input value D(2,1). 1439 !! = 2 the 2nd derivative at TAU(1) is equal to 1440 !! the input value D(2,1). 1441 !! 1442 !! IBCEND : boundary condition indicator at TAU(N). 1443 !! = 0 no boundary condition at TAU(N) is given. 1444 !! In this case, the "not-a-knot condition" 1445 !! is used. 1446 !! = 1 the 1st derivative at TAU(N) is equal to 1447 !! the input value D(2,2). 1448 !! = 2 the 2nd derivative at TAU(N) is equal to 1449 !! the input value D(2,2). 1450 !! 1451 !! D(1,1): value of the function at TAU(1) 1452 !! D(1,1): value of the function at TAU(N) 1453 !! D(2,1): if IBCBEG is 1 (2) it is the value of the 1454 !! 1st (2nd) derivative at TAU(1) 1455 !! D(2,2): if IBCBEG is 1 (2) it is the value of the 1456 !! 1st (2nd) derivative at TAU(N) 1457 !! 1458 !! Output: 1459 !! C(4,N) : contains the polynomial coefficients of the 1460 !! cubic interpolating spline. 1461 !! In the interval interval (TAU(I), TAU(I+1)), 1462 !! the spline F is given by 1463 !! 1464 !! F(X) = 1465 !! C(1,I) + 1466 !! C(2,I) * H + 1467 !! C(3,I) * H^2 / 2 + 1468 !! C(4,I) * H^3 / 6. 1469 !! 1470 !! where H = X - TAU(I). 1471 !! 1472 !!---------------------------------------------------------------------- 1473 INTEGER, INTENT(in ) :: n 1474 REAL(wp), INTENT(in ) :: tau(n) 1475 REAL(wp), INTENT(in ) :: d(2,2) 1476 INTEGER, INTENT(in ) :: ibcbeg, ibcend 1477 REAL(wp) :: c(4,n) 1478 REAL(wp) :: divdf1 1479 REAL(wp) :: divdf3 1480 REAL(wp) :: dtau 1481 REAL(wp) :: g 1482 INTEGER(wp) :: i 1483 !!---------------------------------------------------------------------- 1484 ! Initialise c and copy d values to c 1485 c(:,:) = 0.0D+00 1486 c(1,1) = d(1,1) 1487 c(1,n) = d(1,2) 1488 c(2,1) = d(2,1) 1489 c(2,n) = d(2,2) 1490 ! 1491 ! C(3,*) and C(4,*) are used initially for temporary storage. 1492 ! Store first differences of the TAU sequence in C(3,*). 1493 ! Store first divided difference of data in C(4,*). 1494 DO i = 2, n 1495 c(3,i) = tau(i) - tau(i-1) 1496 c(4,i) = ( c(1,i) - c(1,i-1) ) / c(3,i) 1497 END DO 1498 1499 ! Construct the first equation from the boundary condition 1500 ! at the left endpoint, of the form: 1501 ! 1502 ! C(4,1) * S(1) + C(3,1) * S(2) = C(2,1) 1503 ! 1504 ! IBCBEG = 0: Not-a-knot 1505 IF ( ibcbeg == 0 ) THEN 1506 1507 IF ( n <= 2 ) THEN 1508 c(4,1) = 1._wp 1509 c(3,1) = 1._wp 1510 c(2,1) = 2._wp * c(4,2) 1511 ELSE 1512 c(4,1) = c(3,3) 1513 c(3,1) = c(3,2) + c(3,3) 1514 c(2,1) = ( ( c(3,2) + 2._wp * c(3,1) ) * c(4,2) & 1515 * c(3,3) + c(3,2)**2 * c(4,3) ) / c(3,1) 1516 END IF 1517 1518 ! IBCBEG = 1: derivative specified. 1519 ELSE IF ( ibcbeg == 1 ) then 1520 1521 c(4,1) = 1._wp 1522 c(3,1) = 0._wp 1523 1524 ! IBCBEG = 2: Second derivative prescribed at left end. 1525 ELSE IF ( ibcbeg == 2 ) then 1526 1527 c(4,1) = 2._wp 1528 c(3,1) = 1._wp 1529 c(2,1) = 3._wp * c(4,2) - c(3,2) / 2._wp * c(2,1) 1530 1531 ELSE 1532 WRITE(ctlmes,*) 'CUBSPL - Error, invalid IBCBEG input option!' 1533 CALL ctl_stop( ctlmes ) 1534 END IF 1535 1536 ! If there are interior knots, generate the corresponding 1537 ! equations and carry out the forward pass of Gauss, 1538 ! elimination after which the I-th equation reads: 1539 ! 1540 ! C(4,I) * S(I) + C(3,I) * S(I+1) = C(2,I). 1541 1542 IF ( n > 2 ) THEN 1543 1544 DO i = 2, n-1 1545 g = -c(3,i+1) / c(4,i-1) 1546 c(2,i) = g * c(2,i-1) + 3._wp * & 1547 ( c(3,i) * c(4,i+1) + c(3,i+1) * c(4,i) ) 1548 c(4,i) = g * c(3,i-1) + 2._wp * ( c(3,i) + c(3,i+1)) 1549 END DO 1550 1551 ! Construct the last equation from the second boundary, 1552 ! condition of the form 1553 ! 1554 ! -G * C(4,N-1) * S(N-1) + C(4,N) * S(N) = C(2,N) 1555 ! 1556 ! If 1st der. is prescribed at right end ( ibcend == 1 ), 1557 ! one can go directly to back-substitution, since the C 1558 ! array happens to be set up just right for it at this 1559 ! point. 1560 1561 IF ( ibcend < 1 ) THEN 1562 ! Not-a-knot and 3 <= N, and either 3 < N or also 1563 ! not-a-knot at left end point. 1564 IF ( n /= 3 .OR. ibcbeg /= 0 ) THEN 1565 g = c(3,n-1) + c(3,n) 1566 c(2,n) = ( ( c(3,n) + 2._wp * g ) * c(4,n) * c(3,n-1) & 1567 + c(3,n)**2 * ( c(1,n-1) - c(1,n-2) ) / & 1568 c(3,n-1) ) / g 1569 g = - g / c(4,n-1) 1570 c(4,n) = c(3,n-1) 1571 c(4,n) = c(4,n) + g * c(3,n-1) 1572 c(2,n) = ( g * c(2,n-1) + c(2,n) ) / c(4,n) 1573 ELSE 1574 ! N = 3 and not-a-knot also at left. 1575 c(2,n) = 2._wp * c(4,n) 1576 c(4,n) = 1._wp 1577 g = -1._wp / c(4,n-1) 1578 c(4,n) = c(4,n) - c(3,n-1) / c(4,n-1) 1579 c(2,n) = ( g * c(2,n-1) + c(2,n) ) / c(4,n) 1580 END IF 1581 1582 ELSE IF ( ibcend == 2 ) THEN 1583 ! IBCEND = 2: Second derivative prescribed at right endpoint. 1584 c(2,n) = 3._wp * c(4,n) + c(3,n) / 2._wp * c(2,n) 1585 c(4,n) = 2._wp 1586 g = -1._wp / c(4,n-1) 1587 c(4,n) = c(4,n) - c(3,n-1) / c(4,n-1) 1588 c(2,n) = ( g * c(2,n-1) + c(2,n) ) / c(4,n) 1589 END IF 1590 1591 ELSE 1592 ! N = 2 (assumed to be at least equal to 2!). 1593 1594 IF ( ibcend == 2 ) THEN 1595 1596 c(2,n) = 3._wp * c(4,n) + c(3,n) / 2._wp * c(2,n) 1597 c(4,n) = 2._wp 1598 g = -1._wp / c(4,n-1) 1599 c(4,n) = c(4,n) - c(3,n-1) / c(4,n-1) 1600 c(2,n) = ( g * c(2,n-1) + c(2,n) ) / c(4,n) 1601 1602 ELSE IF ( ibcend == 0 .AND. ibcbeg /= 0 ) THEN 1603 1604 c(2,n) = 2._wp * c(4,n) 1605 c(4,n) = 1._wp 1606 g = -1._wp / c(4,n-1) 1607 c(4,n) = c(4,n) - c(3,n-1) / c(4,n-1) 1608 c(2,n) = ( g * c(2,n-1) + c(2,n) ) / c(4,n) 1609 1610 ELSE IF ( ibcend == 0 .AND. ibcbeg == 0 ) THEN 1611 1612 c(2,n) = c(4,n) 1613 1614 END IF 1615 1616 END IF 1617 1618 ! Back solve the upper triangular system 1619 ! 1620 ! C(4,I) * S(I) + C(3,I) * S(I+1) = B(I) 1621 ! 1622 ! for the slopes C(2,I), given that S(N) is already known. 1623 ! 1624 DO i = n-1, 1, -1 1625 c(2,i) = ( c(2,i) - c(3,i) * c(2,i+1) ) / c(4,i) 1626 END DO 1627 1628 ! Generate cubic coefficients in each interval, that is, the 1629 ! derivatives at its left endpoint, from value and slope at its 1630 ! endpoints. 1631 DO i = 2, n 1632 dtau = c(3,i) 1633 divdf1 = ( c(1,i) - c(1,i-1) ) / dtau 1634 divdf3 = c(2,i-1) + c(2,i) - 2._wp * divdf1 1635 c(3,i-1) = 2._wp * ( divdf1 - c(2,i-1) - divdf3 ) / dtau 1636 c(4,i-1) = 6._wp * divdf3 / dtau**2 1637 END DO 1638 1639 END FUNCTION cub_spl 1640 1641 ! ===================================================================================================== 1642 1643 FUNCTION z_cub_spl(s, tau, c) RESULT ( z_cs ) 1644 !---------------------------------------------------------------------- 1645 ! *** function z_cub_spl *** 1646 ! 1647 ! ** Purpose : evaluate the cubic spline F in the interval 1648 ! (TAU(I), TAU(I+1)). F is given by 1649 ! 1650 ! 1651 ! F(X) = C1 + C2*H + (C3*H^2)/2 + (C4*H^3)/6 1652 ! 1653 ! where H = S-TAU(I). 1654 ! 1655 ! s MUST be positive: 0 <= s <= 1 1656 !--------------------------------------------------------------------- 1657 INTEGER, PARAMETER :: n = 2 1658 REAL(wp), INTENT(in ) :: s 1659 REAL(wp), INTENT(in ) :: tau(n) 1660 REAL(wp), INTENT(in ) :: c(4,n) 1661 REAL(wp) :: z_cs 1662 REAL(wp) :: c1, c2, c3, c4 1663 REAL(wp) :: H 1664 !--------------------------------------------------------------------- 1665 c1 = c(1,1) 1666 c2 = c(2,1) 1667 c3 = c(3,1) 1668 c4 = c(4,1) 1669 1670 H = s - tau(1) 1671 1672 z_cs = c1 + c2 * H + (c3 * H**2._wp)/2._wp + (c4 * H**3._wp)/6._wp 1673 1674 END FUNCTION z_cub_spl 1675 1676 ! ===================================================================================================== 1677 1678 FUNCTION D1z_cub_spl(s, tau, c, kmax) RESULT ( d1z_cs ) 1679 !---------------------------------------------------------------------- 1680 ! *** function D1z_cub_spl *** 1681 ! 1682 ! ** Purpose : evaluate the 1st derivative of cubic spline F 1683 ! in the interval (TAU(I), TAU(I+1)). 1684 ! The 1st derivative D1F is given by 1685 ! 1686 ! 1687 ! D1F(S) = C1 + C2*H + (C3*H^2)/2 + (C4*H^3)/6 1688 ! 1689 ! where H = S-TAU(I). 1690 ! 1691 ! s MUST be positive: 0 <= s <= 1 1692 !--------------------------------------------------------------------- 1693 INTEGER, PARAMETER :: n = 2 1694 REAL(wp), INTENT(in ) :: s 1695 REAL(wp), INTENT(in ) :: tau(n) 1696 REAL(wp), INTENT(in ) :: c(4,n) 1697 INTEGER, INTENT(in ) :: kmax 1698 REAL(wp) :: d1z_cs 1699 REAL(wp) :: c2, c3, c4 1700 REAL(wp) :: H 1701 !--------------------------------------------------------------------- 1702 c2 = c(2,1) / REAL(kmax,wp) 1703 c3 = c(3,1) / REAL(kmax,wp) 1704 c4 = c(4,1) / REAL(kmax,wp) 1705 1706 H = s - tau(1) 1707 1708 d1z_cs = c2 + c3 * H + (c4 * H**2._wp)/2._wp 1709 1710 END FUNCTION D1z_cub_spl 1711 1712 ! ===================================================================================================== 1713 269 END DO 270 END DO 271 ! Computing gdep3w 272 gde3w_0(:,:,1) = 0.5 * e3w_0(:,:,1) 273 DO jk = 2, jpk 274 gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk) 275 END DO 276 ! 277 ! From here equal to sco code - domzgr.F90 line 2079 278 ! Lateral B.C. we apply only for transition zone 279 ! since for s- and zps- areas has already been applied 280 281 CALL lbc_lnk( e3t_0 , 'T', 1._wp ) 282 CALL lbc_lnk( e3u_0 , 'U', 1._wp ) 283 CALL lbc_lnk( e3v_0 , 'V', 1._wp ) 284 CALL lbc_lnk( e3f_0 , 'F', 1._wp ) 285 CALL lbc_lnk( e3w_0 , 'W', 1._wp ) 286 CALL lbc_lnk( e3uw_0, 'U', 1._wp ) 287 CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 288 289 WHERE (e3t_0 (:,:,:).eq.0.0) e3t_0(:,:,:) = 1.0 290 WHERE (e3u_0 (:,:,:).eq.0.0) e3u_0(:,:,:) = 1.0 291 WHERE (e3v_0 (:,:,:).eq.0.0) e3v_0(:,:,:) = 1.0 292 WHERE (e3f_0 (:,:,:).eq.0.0) e3f_0(:,:,:) = 1.0 293 WHERE (e3w_0 (:,:,:).eq.0.0) e3w_0(:,:,:) = 1.0 294 WHERE (e3uw_0 (:,:,:).eq.0.0) e3uw_0(:,:,:) = 1.0 295 WHERE (e3vw_0 (:,:,:).eq.0.0) e3vw_0(:,:,:) = 1.0 296 297 IF( lwp ) WRITE(numout,*) 'Refine mbathy' 298 DO jj = 1, jpj 299 DO ji = 1, jpi 300 IF ( s2z_msk(ji,jj) == 1._wp ) THEN 301 DO jk = 1, jpkm1 302 IF( scobot(ji,jj) >= gdept_0(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) 303 END DO 304 IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0 305 END IF 306 END DO 307 END DO 308 309 END SUBROUTINE zgr_mes_local 1714 310 1715 311 END MODULE zgrmes
Note: See TracChangeset
for help on using the changeset viewer.
