- Timestamp:
- 2018-09-12T17:39:01+02:00 (6 years ago)
- Location:
- NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE
- Files:
-
- 4 edited
-
DIA/diawri.F90 (modified) (1 diff)
-
DOM/domvvl.F90 (modified) (30 diffs)
-
DYN/dynnxt.F90 (modified) (3 diffs)
-
SBC/sbcwave.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/DIA/diawri.F90
r9652 r10116 125 125 ! Output of initial vertical scale factor 126 126 CALL iom_put("e3t_0", e3t_0(:,:,:) ) 127 CALL iom_put("e3u_0", e3 t_0(:,:,:) )128 CALL iom_put("e3v_0", e3 t_0(:,:,:) )127 CALL iom_put("e3u_0", e3u_0(:,:,:) ) 128 CALL iom_put("e3v_0", e3v_0(:,:,:) ) 129 129 ! 130 130 CALL iom_put( "e3t" , e3t_n(:,:,:) ) -
NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/DOM/domvvl.F90
r10060 r10116 8 8 !! 3.3 ! 2011-10 (M. Leclair) totally rewrote domvvl: vvl option includes z_star and z_tilde coordinates 9 9 !! 3.6 ! 2014-11 (P. Mathiot) add ice shelf capability 10 !! 4.0 ! 2018-09 (J. Chanut) improve z_tilde robustness 10 11 !!---------------------------------------------------------------------- 11 12 … … 31 32 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 32 33 USE timing ! Timing 34 USE bdy_oce ! ocean open boundary conditions 35 USE sbcrnf ! river runoff 36 USE dynspg_ts, ONLY: un_adv, vn_adv 33 37 34 38 IMPLICIT NONE … … 44 48 LOGICAL , PUBLIC :: ln_vvl_ztilde = .FALSE. ! ztilde vertical coordinate 45 49 LOGICAL , PUBLIC :: ln_vvl_layer = .FALSE. ! level vertical coordinate 46 LOGICAL , PUBLIC :: ln_vvl_ztilde_as_zstar = .FALSE. ! ztilde vertical coordinate 47 LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor = .FALSE. ! ztilde vertical coordinate 48 LOGICAL , PUBLIC :: ln_vvl_kepe = .FALSE. ! kinetic/potential energy transfer 50 LOGICAL :: ln_vvl_ztilde_as_zstar = .FALSE. ! ztilde vertical coordinate 51 LOGICAL :: ln_vvl_zstar_at_eqtor = .FALSE. ! ztilde vertical coordinate 52 LOGICAL :: ln_vvl_zstar_on_shelf = .FALSE. ! revert to zstar on shelves 53 LOGICAL :: ln_vvl_adv_fct = .FALSE. ! Centred thickness advection 54 LOGICAL :: ln_vvl_adv_cn2 = .TRUE. ! FCT thickness advection 55 LOGICAL :: ln_vvl_dbg = .FALSE. ! debug control prints 56 LOGICAL :: ln_vvl_ramp = .FALSE. ! Ramp on interfaces displacement 57 LOGICAL :: ln_vvl_lap = .FALSE. ! Laplacian thickness diffusion 58 LOGICAL :: ln_vvl_blp = .FALSE. ! Bilaplacian thickness diffusion 59 LOGICAL :: ln_vvl_regrid = .FALSE. ! ensure layer separation 60 LOGICAL :: ll_shorizd = .FALSE. ! Use "shelf horizon depths" 61 LOGICAL :: ln_vvl_kepe = .FALSE. ! kinetic/potential energy transfer 49 62 ! ! conservation: not used yet 50 REAL(wp) :: rn_ahe3 ! thickness diffusion coefficient 51 REAL(wp) :: rn_rst_e3t ! ztilde to zstar restoration timescale [days] 52 REAL(wp) :: rn_lf_cutoff ! cutoff frequency for low-pass filter [days] 53 REAL(wp) :: rn_zdef_max ! maximum fractional e3t deformation 54 LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE. ! debug control prints 63 REAL(wp) :: rn_ahe3_lap ! thickness diffusion coefficient (Laplacian) 64 REAL(wp) :: rn_ahe3_blp ! thickness diffusion coefficient (Bilaplacian) 65 REAL(wp) :: rn_rst_e3t ! ztilde to zstar restoration timescale [days] 66 REAL(wp) :: rn_lf_cutoff ! cutoff frequency for low-pass filter [days] 67 REAL(wp) :: rn_day_ramp ! Duration of linear ramp [days] 68 REAL(wp) :: hsmall=0.01_wp ! small thickness [m] 55 69 56 70 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td ! thickness diffusion transport 57 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf ! low frequency part of hzdivergence71 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_lf, vn_lf, hdivn_lf ! low frequency fluxes and divergence 58 72 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n ! baroclinic scale factors 59 73 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a ! baroclinic scale factors 60 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_e3t ! retoring period for scale factors 61 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_hdv ! retoring period for low freq. divergence 74 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: tildemask ! mask tilde tendency 75 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_e3t ! restoring period for scale factors 76 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_hdv ! restoring period for low freq. divergence 77 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: hsm, dsm ! 78 INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:) :: i_int_bot 62 79 63 80 !! * Substitutions … … 78 95 ALLOCATE( tilde_e3t_b(jpi,jpj,jpk) , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) , & 79 96 & dtilde_e3t_a(jpi,jpj,jpk) , un_td (jpi,jpj,jpk) , vn_td (jpi,jpj,jpk) , & 80 & STAT = dom_vvl_alloc)97 & tildemask(jpi,jpj) , hsm(jpi,jpj) , dsm(jpi,jpj) , i_int_bot(jpi,jpj), STAT = dom_vvl_alloc ) 81 98 IF( lk_mpp ) CALL mpp_sum ( dom_vvl_alloc ) 82 99 IF( dom_vvl_alloc /= 0 ) CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') … … 85 102 ENDIF 86 103 IF( ln_vvl_ztilde ) THEN 87 ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdiv_lf(jpi,jpj,jpk) , STAT= dom_vvl_alloc ) 104 ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdivn_lf(jpi,jpj,jpk) , & 105 & un_lf(jpi,jpj,jpk), vn_lf(jpi,jpj,jpk) , STAT= dom_vvl_alloc ) 88 106 IF( lk_mpp ) CALL mpp_sum ( dom_vvl_alloc ) 89 107 IF( dom_vvl_alloc /= 0 ) CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') … … 117 135 INTEGER :: ji, jj, jk 118 136 INTEGER :: ii0, ii1, ij0, ij1 119 REAL(wp):: zcoef 137 REAL(wp):: zcoef, zwgt, ztmp, zhmin, zhmax 120 138 !!---------------------------------------------------------------------- 121 139 ! … … 129 147 IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) 130 148 ! 131 ! ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv _lf149 ! ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdivn_lf 132 150 CALL dom_vvl_rst( nit000, 'READ' ) 133 151 e3t_a(:,:,jpk) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all … … 199 217 200 218 ! !== z_tilde coordinate case ==! (Restoring frequencies) 219 tildemask(:,:) = 1._wp 220 201 221 IF( ln_vvl_ztilde ) THEN 202 222 !!gm : idea: add here a READ in a file of custumized restoring frequency 203 ! ! Values in days provided via the namelist 204 ! ! use rsmall to avoid possible division by zero errors with faulty settings 205 frq_rst_e3t(:,:) = 2._wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.0_wp ) 206 frq_rst_hdv(:,:) = 2._wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp ) 207 ! 208 IF( ln_vvl_ztilde_as_zstar ) THEN ! z-star emulation using z-tile 209 frq_rst_e3t(:,:) = 0._wp !Ignore namelist settings 210 frq_rst_hdv(:,:) = 1._wp / rdt 223 ! Values in days provided via the namelist; use rsmall to avoid possible division by zero errors with faulty settings 224 frq_rst_e3t(:,:) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.0_wp ) 225 frq_rst_hdv(:,:) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp ) 226 ! 227 IF( ln_vvl_ztilde_as_zstar ) THEN 228 ! Ignore namelist settings and use these next two to emulate z-star using z-tilde 229 frq_rst_e3t(:,:) = 0.0_wp 230 frq_rst_hdv(:,:) = 1.0_wp / rdt 231 tildemask(:,:) = 0._wp 211 232 ENDIF 212 IF ( ln_vvl_zstar_at_eqtor ) THEN ! use z-star in vicinity of the Equator 233 234 IF ( ln_vvl_zstar_at_eqtor ) THEN 213 235 DO jj = 1, jpj 214 236 DO ji = 1, jpi 215 !!gm case |gphi| >= 6 degrees is useless initialized just above by default 237 !!gm case |gphi| >= 6 degrees is useless initialized just above by default 216 238 IF( ABS(gphit(ji,jj)) >= 6.) THEN 217 239 ! values outside the equatorial band and transition zone (ztilde) 218 240 frq_rst_e3t(ji,jj) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.e0_wp ) 219 frq_rst_hdv(ji,jj) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 220 ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN ! Equator strip ==> z-star 241 ! frq_rst_hdv(ji,jj) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 242 243 ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN 221 244 ! values inside the equatorial band (ztilde as zstar) 222 245 frq_rst_e3t(ji,jj) = 0.0_wp 223 frq_rst_hdv(ji,jj) = 1.0_wp / rdt 224 ELSE ! transition band (2.5 to 6 degrees N/S) 225 ! ! (linearly transition from z-tilde to z-star) 246 ! frq_rst_hdv(ji,jj) = 1.0_wp / rdt 247 tildemask(ji,jj) = 0._wp 248 ELSE 249 ! values in the transition band (linearly vary from ztilde to ztilde as zstar values) 226 250 frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp & 227 251 & * ( 1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & 228 252 & * 180._wp / 3.5_wp ) ) 229 frq_rst_hdv(ji,jj) = (1.0_wp / rdt) & 230 & + ( frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp & 231 & * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & 232 & * 180._wp / 3.5_wp ) ) 253 ! frq_rst_hdv(ji,jj) = (1.0_wp / rdt) & 254 ! & + ( frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp & 255 ! & * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & 256 ! & * 180._wp / 3.5_wp ) ) 257 tildemask(ji,jj) = 0.5_wp * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & 258 & * 180._wp / 3.5_wp ) ) 233 259 ENDIF 234 260 END DO 235 261 END DO 236 IF( cn_cfg == "orca" .AND. nn_cfg == 3 ) THEN ! ORCA2: Suppress ztilde in the Foxe Basin for ORCA2237 ii0 = 103 ; ii1 = 111238 ij0 = 128 ; ij1 = 135 ;239 frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.0_wp240 frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0_wp / rdt241 ENDIF242 262 ENDIF 263 ! 264 IF ( ln_vvl_zstar_on_shelf ) THEN 265 zhmin = 50._wp 266 zhmax = 100._wp 267 DO jj = 1, jpj 268 DO ji = 1, jpi 269 zwgt = 1._wp 270 IF(( ht_0(ji,jj)>zhmin).AND.(ht_0(ji,jj) <=zhmax)) THEN 271 zwgt = (ht_0(ji,jj)-zhmin)/(zhmax-zhmin) 272 ELSEIF ( ht_0(ji,jj)<=zhmin) THEN 273 zwgt = 0._wp 274 ENDIF 275 frq_rst_e3t(ji,jj) = MIN(frq_rst_e3t(ji,jj), frq_rst_e3t(ji,jj)*zwgt) 276 tildemask(ji,jj) = MIN(tildemask(ji,jj), zwgt) 277 END DO 278 END DO 279 ENDIF 280 ! 281 ztmp = MAXVAL( frq_rst_hdv(:,:) ) 282 IF( lk_mpp ) CALL mpp_max( ztmp ) ! max over the global domain 283 ! 284 IF ( (ztmp*rdt) > 1._wp) CALL ctl_stop( 'dom_vvl_init: rn_lf_cuttoff is too small' ) 285 ! 243 286 ENDIF 244 287 ! … … 256 299 IF( ln_vvl_ztilde ) THEN ! z_tilde case ! 257 300 ! ! ------------ ! 258 CALL iom_set_rstw_var_active('hdiv_lf') 301 CALL iom_set_rstw_var_active('un_lf') 302 CALL iom_set_rstw_var_active('vn_lf') 303 CALL iom_set_rstw_var_active('hdivn_lf') 259 304 ENDIF 260 305 ! … … 279 324 !! to the "baroclinic" level thickness. 280 325 !! 281 !! ** Action : - hdiv _lf: restoring towards full baroclinic divergence in z_tilde case326 !! ** Action : - hdivn_lf : restoring towards full baroclinic divergence in z_tilde case 282 327 !! - tilde_e3t_a: after increment of vertical scale factor 283 328 !! in z_tilde case … … 289 334 INTEGER, INTENT( in ), OPTIONAL :: kcall ! optional argument indicating call sequence 290 335 ! 336 LOGICAL :: ll_do_bclinic ! local logical 291 337 INTEGER :: ji, jj, jk ! dummy loop indices 338 INTEGER :: ib, ib_bdy, ip, jp ! " " " 292 339 INTEGER , DIMENSION(3) :: ijk_max, ijk_min ! temporary integers 293 REAL(wp) :: z2dt, z_tmin, z_tmax ! local scalars 294 LOGICAL :: ll_do_bclinic ! local logical 340 INTEGER :: ncall 341 REAL(wp) :: z2dt , z_tmin, z_tmax! local scalars 342 REAL(wp) :: zalpha, zwgt ! temporary scalars 343 REAL(wp) :: zdu, zdv, zramp 344 REAL(wp) :: ztra, zbtr, ztout, ztin, zfac, zmsku, zmskv 295 345 REAL(wp), DIMENSION(jpi,jpj) :: zht, z_scale, zwu, zwv, zhdiv 296 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t 346 REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t, ztu, ztv 297 347 !!---------------------------------------------------------------------- 298 348 ! … … 308 358 309 359 ll_do_bclinic = .TRUE. 360 ncall = 1 361 310 362 IF( PRESENT(kcall) ) THEN 311 IF( kcall == 2 .AND. ln_vvl_ztilde ) ll_do_bclinic = .FALSE. 363 ! comment line below if tilda coordinate has to be computed at each call 364 ! This is not working yet 365 ! IF (kcall == 2 .AND. (ln_vvl_ztilde.OR.ln_vvl_layer) ) ll_do_bclinic = .FALSE. 366 ncall = kcall 367 ENDIF 368 369 IF( neuler == 0 .AND. kt == nit000 ) THEN 370 z2dt = rdt 371 ELSE 372 z2dt = 2.0_wp * rdt 312 373 ENDIF 313 374 314 375 ! ******************************* ! 315 ! After acale factors at t-points !376 ! After scale factors at t-points ! 316 377 ! ******************************* ! 317 378 ! ! --------------------------------------------- ! … … 325 386 END DO 326 387 ! 327 IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate ! 328 ! ! ------baroclinic part------ ! 329 ! I - initialization 330 ! ================== 331 332 ! 1 - barotropic divergence 333 ! ------------------------- 334 zhdiv(:,:) = 0._wp 335 zht(:,:) = 0._wp 388 IF( (ln_vvl_ztilde.OR.ln_vvl_layer) .AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate ! 389 ! ! ------baroclinic part------ ! 390 tilde_e3t_a(:,:,:) = 0.0_wp ! tilde_e3t_a used to store tendency terms 391 un_td(:,:,:) = 0.0_wp ! Transport corrections 392 vn_td(:,:,:) = 0.0_wp 393 394 zhdiv(:,:) = 0. 336 395 DO jk = 1, jpkm1 337 396 zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk) 338 zht (:,:) = zht (:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 339 END DO 340 zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) ) 341 342 ! 2 - Low frequency baroclinic horizontal divergence (z-tilde case only) 343 ! -------------------------------------------------- 344 IF( ln_vvl_ztilde ) THEN 345 IF( kt > nit000 ) THEN 397 END DO 398 zhdiv(:,:) = zhdiv(:,:) / ( ht_0(:,:) + sshn(:,:) + 1._wp - tmask_i(:,:) ) 399 400 ze3t(:,:,:) = 0._wp 401 IF( ln_rnf ) THEN 402 CALL sbc_rnf_div( ze3t ) ! runoffs 403 DO jk=1,jpkm1 404 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ze3t(:,:,jk) * e3t_n(:,:,jk) 405 END DO 406 ENDIF 407 408 ! Thickness advection: 409 ! -------------------- 410 ! Set advection velocities and source term 411 IF ( ln_vvl_ztilde ) THEN 412 IF ( ncall==1 ) THEN 413 zalpha = rdt * 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp ) 346 414 DO jk = 1, jpkm1 347 hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:) & 348 & * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) ) 349 END DO 415 ztu(:,:,jk) = un(:,:,jk) * e3u_n(:,:,jk) * e2u(:,:) 416 ztv(:,:,jk) = vn(:,:,jk) * e3v_n(:,:,jk) * e1v(:,:) 417 ze3t(:,:,jk) = -tilde_e3t_a(:,:,jk) - e3t_n(:,:,jk) * zhdiv(:,:) 418 END DO 419 ! 420 un_lf(:,:,:) = un_lf(:,:,:) * (1._wp - zalpha) + zalpha * ztu(:,:,:) 421 vn_lf(:,:,:) = vn_lf(:,:,:) * (1._wp - zalpha) + zalpha * ztv(:,:,:) 422 hdivn_lf(:,:,:) = hdivn_lf(:,:,:) * (1._wp - zalpha) + zalpha * ze3t(:,:,:) 423 ! 2nd order filtering: 424 ! un_lf2(:,:,:) = un_lf2(:,:,:) * (1._wp - zalpha) + zalpha * ztu(:,:,:) 425 ! vn_lf2(:,:,:) = vn_lf2(:,:,:) * (1._wp - zalpha) + zalpha * ztv(:,:,:) 426 ! hdivn_lf2(:,:,:) = hdivn_lf2(:,:,:) * (1._wp - zalpha) + zalpha * ze3t(:,:,:) 427 ! un_lf(:,:,:) = un_lf(:,:,:) * (1._wp - zalpha) + zalpha * un_lf2(:,:,:) 428 ! vn_lf(:,:,:) = vn_lf(:,:,:) * (1._wp - zalpha) + zalpha * vn_lf2(:,:,:) 429 ! hdivn_lf(:,:,:) = hdivn_lf(:,:,:) * (1._wp - zalpha) + zalpha * hdivn_lf2(:,:,:) 350 430 ENDIF 431 ! 432 DO jk = 1, jpkm1 433 ztu(:,:,jk) = (un(:,:,jk)-un_lf(:,:,jk)/e3u_n(:,:,jk)*r1_e2u(:,:))*umask(:,:,jk) 434 ztv(:,:,jk) = (vn(:,:,jk)-vn_lf(:,:,jk)/e3v_n(:,:,jk)*r1_e1v(:,:))*vmask(:,:,jk) 435 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) + hdivn_lf(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk) 436 END DO 437 438 ! 439 ELSEIF ( ln_vvl_layer ) THEN 440 ! 441 DO jk = 1, jpkm1 442 ztu(:,:,jk) = un(:,:,jk) 443 ztv(:,:,jk) = vn(:,:,jk) 444 END DO 445 ! 351 446 ENDIF 352 353 ! II - after z_tilde increments of vertical scale factors 354 ! ======================================================= 355 tilde_e3t_a(:,:,:) = 0._wp ! tilde_e3t_a used to store tendency terms 356 357 ! 1 - High frequency divergence term 358 ! ---------------------------------- 359 IF( ln_vvl_ztilde ) THEN ! z_tilde case 447 ! 448 ! Block fluxes through small layers: 449 DO jk=1,jpkm1 450 DO ji = 1, jpi 451 DO jj= 1, jpj 452 zmsku = 0.5_wp * ( 1._wp + SIGN(1._wp, e3u_n(ji,jj,jk) - hsmall) ) 453 un_td(ji,jj,jk) = un_td(ji,jj,jk) - (1. - zmsku) * un(ji,jj,jk) * e3u_n(ji,jj,jk) * e2u(ji,jj) 454 ztu(ji,jj,jk) = zmsku * ztu(ji,jj,jk) 455 IF ( ln_vvl_ztilde ) un_lf(ji,jj,jk) = zmsku * un_lf(ji,jj,jk) 456 ! 457 zmskv = 0.5_wp * ( 1._wp + SIGN(1._wp, e3v_n(ji,jj,jk) - hsmall) ) 458 vn_td(ji,jj,jk) = vn_td(ji,jj,jk) - (1. - zmskv) * vn(ji,jj,jk) * e3v_n(ji,jj,jk) * e1v(ji,jj) 459 ztv(ji,jj,jk) = zmskv * ztv(ji,jj,jk) 460 IF ( ln_vvl_ztilde ) vn_lf(ji,jj,jk) = zmskv * vn_lf(ji,jj,jk) 461 END DO 462 END DO 463 END DO 464 ! 465 ! Do advection 466 IF (ln_vvl_adv_fct) THEN 467 CALL dom_vvl_adv_fct( kt, tilde_e3t_a, ztu, ztv ) 468 ! 469 ELSEIF (ln_vvl_adv_cn2) THEN 360 470 DO jk = 1, jpkm1 361 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 362 END DO 363 ELSE ! layer case 471 DO jj = 2, jpjm1 472 DO ji = fs_2, fs_jpim1 ! vector opt. 473 tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) & 474 & -( e2u(ji,jj)*e3u_n(ji,jj,jk) * ztu(ji,jj,jk) - e2u(ji-1,jj )*e3u_n(ji-1,jj ,jk) * ztu(ji-1,jj ,jk) & 475 & + e1v(ji,jj)*e3v_n(ji,jj,jk) * ztv(ji,jj,jk) - e1v(ji ,jj-1)*e3v_n(ji ,jj-1,jk) * ztv(ji ,jj-1,jk) ) & 476 & / ( e1t(ji,jj) * e2t(ji,jj) ) 477 END DO 478 END DO 479 END DO 480 ENDIF 481 ! 482 ! Thickness anomaly diffusion: 483 ! ----------------------------- 484 ztu(:,:,:) = 0.0_wp 485 ztv(:,:,:) = 0.0_wp 486 487 ze3t(:,:,1) = 0.e0 488 DO jj = 1, jpj 489 DO ji = 1, jpi 490 DO jk = 2, jpk 491 ze3t(ji,jj,jk) = ze3t(ji,jj,jk-1) + tilde_e3t_b(ji,jj,jk-1) * tmask(ji,jj,jk-1) 492 END DO 493 END DO 494 END DO 495 496 IF ( ln_vvl_blp ) THEN ! Bilaplacian 364 497 DO jk = 1, jpkm1 365 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 498 DO jj = 1, jpjm1 ! First derivative (gradient) 499 DO ji = 1, fs_jpim1 ! vector opt. 500 ! ztu(ji,jj,jk) = umask(ji,jj,jk) * e2_e1u(ji,jj) & 501 ! & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj ,jk) ) 502 ! ztv(ji,jj,jk) = vmask(ji,jj,jk) * e1_e2v(ji,jj) & 503 ! & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji ,jj+1,jk) ) 504 ztu(ji,jj,jk) = umask(ji,jj,jk) * e2_e1u(ji,jj) & 505 & * ( ze3t(ji,jj,jk) - ze3t(ji+1,jj ,jk) ) 506 ztv(ji,jj,jk) = vmask(ji,jj,jk) * e1_e2v(ji,jj) & 507 & * ( ze3t(ji,jj,jk) - ze3t(ji ,jj+1,jk) ) 508 END DO 509 END DO 510 511 DO jj = 2, jpjm1 ! Second derivative (divergence) time the eddy diffusivity coefficient 512 DO ji = fs_2, fs_jpim1 ! vector opt. 513 zht(ji,jj) = rn_ahe3_blp * r1_e1e2t(ji,jj) * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) & 514 & + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) 515 END DO 516 END DO 517 518 ! Open boundary conditions: 519 IF ( ln_bdy ) THEN 520 DO ib_bdy=1, nb_bdy 521 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(1) 522 ji = idx_bdy(ib_bdy)%nbi(ib,1) 523 jj = idx_bdy(ib_bdy)%nbj(ib,1) 524 zht(ji,jj) = 0._wp 525 END DO 526 END DO 527 END IF 528 529 CALL lbc_lnk( zht, 'T', 1. ) ! Lateral boundary conditions (unchanged sgn) 530 531 DO jj = 1, jpjm1 ! third derivative (gradient) 532 DO ji = 1, fs_jpim1 ! vector opt. 533 ztu(ji,jj,jk) = umask(ji,jj,jk) * e2_e1u(ji,jj) * ( zht(ji+1,jj ) - zht(ji,jj) ) 534 ztv(ji,jj,jk) = vmask(ji,jj,jk) * e1_e2v(ji,jj) * ( zht(ji ,jj+1) - zht(ji,jj) ) 535 END DO 536 END DO 366 537 END DO 367 538 ENDIF 368 539 369 ! 2 - Restoring term (z-tilde case only) 370 ! ------------------ 371 IF( ln_vvl_ztilde ) THEN 372 DO jk = 1, jpk 373 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk) 374 END DO 375 ENDIF 376 377 ! 3 - Thickness diffusion term 378 ! ---------------------------- 379 zwu(:,:) = 0._wp 380 zwv(:,:) = 0._wp 381 DO jk = 1, jpkm1 ! a - first derivative: diffusive fluxes 382 DO jj = 1, jpjm1 383 DO ji = 1, fs_jpim1 ! vector opt. 384 un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * e2_e1u(ji,jj) & 385 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj ,jk) ) 386 vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * e1_e2v(ji,jj) & 387 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji ,jj+1,jk) ) 388 zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 389 zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 390 END DO 391 END DO 392 END DO 393 DO jj = 1, jpj ! b - correction for last oceanic u-v points 394 DO ji = 1, jpi 395 un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 396 vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 397 END DO 398 END DO 399 DO jk = 1, jpkm1 ! c - second derivative: divergence of diffusive fluxes 540 IF ( ln_vvl_lap ) THEN ! Laplacian 541 DO jk = 1, jpkm1 ! First derivative (gradient) 542 DO jj = 1, jpjm1 543 DO ji = 1, fs_jpim1 ! vector opt. 544 ! zdu = rn_ahe3_lap * umask(ji,jj,jk) * e2_e1u(ji,jj) & 545 ! & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj ,jk) ) 546 ! zdv = rn_ahe3_lap * vmask(ji,jj,jk) * e1_e2v(ji,jj) & 547 ! & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji ,jj+1,jk) ) 548 zdu = rn_ahe3_lap * umask(ji,jj,jk) * e2_e1u(ji,jj) & 549 & * ( ze3t(ji,jj,jk) - ze3t(ji+1,jj ,jk) ) 550 zdv = rn_ahe3_lap * vmask(ji,jj,jk) * e1_e2v(ji,jj) & 551 & * ( ze3t(ji,jj,jk) - ze3t(ji ,jj+1,jk) ) 552 ztu(ji,jj,jk) = ztu(ji,jj,jk) + zdu 553 ztv(ji,jj,jk) = ztv(ji,jj,jk) + zdv 554 END DO 555 END DO 556 END DO 557 ENDIF 558 559 ! divergence of diffusive fluxes 560 DO jk = 1, jpkm1 400 561 DO jj = 2, jpjm1 401 562 DO ji = fs_2, fs_jpim1 ! vector opt. 402 tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + ( un_td(ji-1,jj ,jk) - un_td(ji,jj,jk) & 403 & + vn_td(ji ,jj-1,jk) - vn_td(ji,jj,jk) & 563 ! tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + ( ztu(ji-1,jj ,jk) - ztu(ji,jj,jk) & 564 ! & + ztv(ji ,jj-1,jk) - ztv(ji,jj,jk) & 565 ! & ) * r1_e1e2t(ji,jj) 566 un_td(ji,jj,jk) = un_td(ji,jj,jk) + ztu(ji,jj,jk+1) - ztu(ji,jj,jk ) 567 vn_td(ji,jj,jk) = vn_td(ji,jj,jk) + ztv(ji,jj,jk+1) - ztv(ji,jj,jk ) 568 tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + ( ztu(ji-1,jj ,jk+1) - ztu(ji,jj,jk+1) & 569 & +ztv(ji ,jj-1,jk+1) - ztv(ji,jj,jk+1) & 570 & -ztu(ji-1,jj ,jk ) + ztu(ji,jj,jk ) & 571 & -ztv(ji ,jj-1,jk ) + ztv(ji,jj,jk ) & 404 572 & ) * r1_e1e2t(ji,jj) 405 573 END DO 406 574 END DO 407 575 END DO 408 ! ! d - thickness diffusion transport: boundary conditions 409 ! (stored for tracer advction and continuity equation) 410 CALL lbc_lnk_multi( un_td , 'U' , -1._wp, vn_td , 'V' , -1._wp) 411 412 ! 4 - Time stepping of baroclinic scale factors 413 ! --------------------------------------------- 576 577 578 ! un_td(:,:,:) = un_td(:,:,:) + ztu(:,:,:) 579 ! vn_td(:,:,:) = vn_td(:,:,:) + ztv(:,:,:) 580 581 CALL lbc_lnk( un_td , 'U' , -1.) 582 CALL lbc_lnk( vn_td , 'V' , -1.) 583 ! 584 CALL dom_vvl_ups_cor( kt, tilde_e3t_a, un_td, vn_td ) 585 586 ! IF ( ln_vvl_ztilde ) THEN 587 ! ztu(:,:,:) = -un_lf(:,:,:) 588 ! ztv(:,:,:) = -vn_lf(:,:,:) 589 ! CALL dom_vvl_ups_cor( kt, tilde_e3t_a, ztu, ztv ) 590 ! ENDIF 591 ! 592 ! Remove "external thickness" tendency: 593 DO jk = 1, jpkm1 594 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) + e3t_n(:,:,jk) * zhdiv(:,:) 595 END DO 596 414 597 ! Leapfrog time stepping 415 598 ! ~~~~~~~~~~~~~~~~~~~~~~ 416 IF( neuler == 0 .AND. kt == nit000 ) THEN 417 z2dt = rdt 418 ELSE 419 z2dt = 2.0_wp * rdt 599 zramp = 1._wp 600 IF ((.NOT.ln_rstart).AND.ln_vvl_ramp) zramp = MIN(MAX( ((kt-nit000)*rdt)/(rn_day_ramp*rday),0._wp),1._wp) 601 602 DO jk=1,jpkm1 603 tilde_e3t_a(:,:,jk) = tilde_e3t_b(:,:,jk) + z2dt * tmask(:,:,jk) * tilde_e3t_a(:,:,jk) & 604 & * tildemask(:,:) * zramp 605 END DO 606 607 ! Ensure layer separation: 608 ! ~~~~~~~~~~~~~~~~~~~~~~~~ 609 IF ( ln_vvl_regrid ) CALL dom_vvl_regrid( kt ) 610 611 ! Boundary conditions: 612 ! ~~~~~~~~~~~~~~~~~~~~ 613 IF ( ln_bdy ) THEN 614 DO ib_bdy = 1, nb_bdy 615 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(1) 616 !! DO ib = 1, idx_bdy(ib_bdy)%nblen(1) 617 ji = idx_bdy(ib_bdy)%nbi(ib,1) 618 jj = idx_bdy(ib_bdy)%nbj(ib,1) 619 zwgt = idx_bdy(ib_bdy)%nbw(ib,1) 620 ip = bdytmask(ji+1,jj ) - bdytmask(ji-1,jj ) 621 jp = bdytmask(ji ,jj+1) - bdytmask(ji ,jj-1) 622 DO jk = 1, jpkm1 623 tilde_e3t_a(ji,jj,jk) = 0.e0 624 !! tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) * (1._wp - zwgt) 625 !! tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji+ip,jj+jp,jk) * tmask(ji+ip,jj+jp,jk) 626 END DO 627 END DO 628 END DO 420 629 ENDIF 421 CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1._wp ) 422 tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:) 423 424 ! Maximum deformation control 425 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 426 ze3t(:,:,jpk) = 0._wp 630 631 CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1. ) 632 633 ENDIF 634 635 IF( ln_vvl_ztilde.AND.( ncall==1)) THEN 636 zalpha = rdt * 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp ) 637 ! 638 ! divergence of diffusive fluxes 427 639 DO jk = 1, jpkm1 428 ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) 429 END DO 430 z_tmax = MAXVAL( ze3t(:,:,:) ) 431 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 432 z_tmin = MINVAL( ze3t(:,:,:) ) 640 DO jj = 2, jpjm1 641 DO ji = fs_2, fs_jpim1 ! vector opt. 642 ze3t(ji,jj,jk) = ( un_td(ji,jj,jk) - un_td(ji-1,jj ,jk) & 643 & + vn_td(ji,jj,jk) - vn_td(ji ,jj-1,jk) & 644 & ) / ( e1t(ji,jj) * e2t(ji,jj) ) 645 END DO 646 END DO 647 END DO 648 CALL lbc_lnk( ze3t(:,:,:), 'T', 1. ) 649 hdivn_lf(:,:,:) = hdivn_lf(:,:,:) + zalpha * ze3t(:,:,:) 650 ! 2nd order filtering: 651 ! hdivn_lf2(:,:,:) = hdivn_lf2(:,:,:) + zalpha * ze3t(:,:,:) 652 ! hdivn_lf(:,:,:) = hdivn_lf(:,:,:) + zalpha * zalpha * ze3t(:,:,:) 653 ENDIF 654 655 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde or layer coordinate ! 656 ! ! ---baroclinic part--------- ! 657 DO jk = 1, jpkm1 658 e3t_a(:,:,jk) = e3t_a(:,:,jk) + (tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)) 659 END DO 660 ENDIF 661 662 IF( ln_vvl_dbg.AND.(ln_vvl_ztilde .OR. ln_vvl_layer) ) THEN ! - ML - test: control prints for debuging 663 ! 664 zht(:,:) = 0.0_wp 665 DO jk = 1, jpkm1 666 zht(:,:) = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk) 667 END DO 668 IF( lwp ) WRITE(numout, *) 'kt = ', kt 669 IF( lwp ) WRITE(numout, *) 'ncall = ', ncall 670 IF( lwp ) WRITE(numout, *) 'll_do_bclinic', ll_do_bclinic 671 IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 672 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ), mask = tmask(:,:,1) == 1.e0 ) 673 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 674 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax 675 END IF 676 ! 677 z_tmin = MINVAL( e3t_n(:,:,:) ) 433 678 IF( lk_mpp ) CALL mpp_min( z_tmin ) ! min over the global domain 434 ! - ML - test: for the moment, stop simulation for too large e3_t variations435 IF ( ( z_tmax > rn_zdef_max ) .OR. ( z_tmin < - rn_zdef_max )) THEN679 IF( lwp ) WRITE(numout, *) kt,' MINVAL(e3t_n) =', z_tmin 680 IF ( z_tmin .LE. 0._wp ) THEN 436 681 IF( lk_mpp ) THEN 437 CALL mpp_maxloc( ze3t, tmask, z_tmax, ijk_max(1), ijk_max(2), ijk_max(3) ) 438 CALL mpp_minloc( ze3t, tmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) ) 682 CALL mpp_minloc(e3t_n(:,:,:), tmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) ) 439 683 ELSE 440 ijk_max = MAXLOC( ze3t(:,:,:) ) 441 ijk_max(1) = ijk_max(1) + nimpp - 1 442 ijk_max(2) = ijk_max(2) + njmpp - 1 443 ijk_min = MINLOC( ze3t(:,:,:) ) 684 ijk_min = MINLOC( e3t_n(:,:,:) ) 444 685 ijk_min(1) = ijk_min(1) + nimpp - 1 445 686 ijk_min(2) = ijk_min(2) + njmpp - 1 446 687 ENDIF 447 688 IF (lwp) THEN 448 WRITE(numout, *) 'MAX( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax 449 WRITE(numout, *) 'at i, j, k=', ijk_max 450 WRITE(numout, *) 'MIN( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin 689 WRITE(numout, *) 'Negative scale factor, e3t_n =', z_tmin 451 690 WRITE(numout, *) 'at i, j, k=', ijk_min 452 CALL ctl_warn('MAX( ABS( tilde_e3t_a(:,:,:) ) / e3t_0(:,:,:) ) too high') 691 CALL ctl_stop('dom_vvl_sf_nxt: Negative scale factor') 692 ENDIF 693 ENDIF 694 ! 695 z_tmin = MINVAL( e3u_n(:,:,:)) 696 IF( lk_mpp ) CALL mpp_min( z_tmin ) ! min over the global domain 697 IF( lwp ) WRITE(numout, *) kt,' MINVAL(e3u_n) =', z_tmin 698 IF ( z_tmin .LE. 0._wp ) THEN 699 IF( lk_mpp ) THEN 700 CALL mpp_minloc(e3u_n(:,:,:), umask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) ) 701 ELSE 702 ijk_min = MINLOC( e3u_n(:,:,:) ) 703 ijk_min(1) = ijk_min(1) + nimpp - 1 704 ijk_min(2) = ijk_min(2) + njmpp - 1 705 ENDIF 706 IF (lwp) THEN 707 WRITE(numout, *) 'Negative scale factor, e3u_n =', z_tmin 708 WRITE(numout, *) 'at i, j, k=', ijk_min 709 CALL ctl_stop('dom_vvl_sf_nxt: Negative scale factor') 710 ENDIF 711 ENDIF 712 ! 713 z_tmin = MINVAL( e3v_n(:,:,:) ) 714 IF( lk_mpp ) CALL mpp_min( z_tmin ) ! min over the global domain 715 IF( lwp ) WRITE(numout, *) kt,' MINVAL(e3v_n) =', z_tmin 716 IF ( z_tmin .LE. 0._wp ) THEN 717 IF( lk_mpp ) THEN 718 CALL mpp_minloc(e3v_n(:,:,:), vmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) ) 719 ELSE 720 ijk_min = MINLOC( e3v_n(:,:,:), mask = vmask(:,:,:) == 1.e0 ) 721 ijk_min(1) = ijk_min(1) + nimpp - 1 722 ijk_min(2) = ijk_min(2) + njmpp - 1 723 ENDIF 724 IF (lwp) THEN 725 WRITE(numout, *) 'Negative scale factor, e3v_n =', z_tmin 726 WRITE(numout, *) 'at i, j, k=', ijk_min 727 CALL ctl_stop('dom_vvl_sf_nxt: Negative scale factor') 728 ENDIF 729 ENDIF 730 ! 731 z_tmin = MINVAL( e3f_n(:,:,:)) 732 IF( lk_mpp ) CALL mpp_min( z_tmin ) ! min over the global domain 733 IF( lwp ) WRITE(numout, *) kt,' MINVAL(e3f_n) =', z_tmin 734 IF ( z_tmin .LE. 0._wp ) THEN 735 IF( lk_mpp ) THEN 736 CALL mpp_minloc(e3f_n(:,:,:), fmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) ) 737 ELSE 738 ijk_min = MINLOC( e3f_n(:,:,:) ) 739 ijk_min(1) = ijk_min(1) + nimpp - 1 740 ijk_min(2) = ijk_min(2) + njmpp - 1 741 ENDIF 742 IF (lwp) THEN 743 WRITE(numout, *) 'Negative scale factor, e3f_n =', z_tmin 744 WRITE(numout, *) 'at i, j, k=', ijk_min 745 CALL ctl_stop('dom_vvl_sf_nxt: Negative scale factor') 453 746 ENDIF 454 747 ENDIF 455 ! - ML - end test456 ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below457 tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:), rn_zdef_max * e3t_0(:,:,:) )458 tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) )459 460 !461 ! "tilda" change in the after scale factor462 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~463 DO jk = 1, jpkm1464 dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)465 END DO466 ! III - Barotropic repartition of the sea surface height over the baroclinic profile467 ! ==================================================================================468 ! add ( ssh increment + "baroclinicity error" ) proportionly to e3t(n)469 ! - ML - baroclinicity error should be better treated in the future470 ! i.e. locally and not spread over the water column.471 ! (keep in mind that the idea is to reduce Eulerian velocity as much as possible)472 zht(:,:) = 0.473 DO jk = 1, jpkm1474 zht(:,:) = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk)475 END DO476 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )477 DO jk = 1, jpkm1478 dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk)479 END DO480 481 ENDIF482 483 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde or layer coordinate !484 ! ! ---baroclinic part--------- !485 DO jk = 1, jpkm1486 e3t_a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk)487 END DO488 ENDIF489 490 IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN ! - ML - test: control prints for debuging491 !492 IF( lwp ) WRITE(numout, *) 'kt =', kt493 IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN494 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) )495 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain496 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax497 END IF498 748 ! 499 749 zht(:,:) = 0.0_wp … … 503 753 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) ) 504 754 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 505 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn -SUM(e3t_n))) =', z_tmax755 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn -SUM(e3t_n))) =', z_tmax 506 756 ! 507 757 zht(:,:) = 0.0_wp 508 758 DO jk = 1, jpkm1 509 zht(:,:) = zht(:,:) + e3t_a(:,:,jk) * tmask(:,:,jk) 510 END DO 511 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) ) 759 zht(:,:) = zht(:,:) + e3u_n(:,:,jk) * umask(:,:,jk) 760 END DO 761 zwu(:,:) = 0._wp 762 DO jj = 1, jpjm1 763 DO ji = 1, fs_jpim1 ! vector opt. 764 zwu(ji,jj) = 0.5_wp * umask(ji,jj,1) * r1_e1e2u(ji,jj) & 765 & * ( e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji+1,jj) * sshn(ji+1,jj) ) 766 END DO 767 END DO 768 CALL lbc_lnk( zwu(:,:), 'U', 1._wp ) 769 z_tmax = MAXVAL( ssumask(:,:) * ssumask(:,:) * ABS( hu_0(:,:) + zwu(:,:) - zht(:,:) ) ) 512 770 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 513 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(h t_0+ssha-SUM(e3t_a))) =', z_tmax771 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(hu_0+sshu_n-SUM(e3u_n))) =', z_tmax 514 772 ! 515 773 zht(:,:) = 0.0_wp 516 774 DO jk = 1, jpkm1 517 zht(:,:) = zht(:,:) + e3t_b(:,:,jk) * tmask(:,:,jk) 518 END DO 519 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) ) 775 zht(:,:) = zht(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk) 776 END DO 777 zwv(:,:) = 0._wp 778 DO jj = 1, jpjm1 779 DO ji = 1, fs_jpim1 ! vector opt. 780 zwv(ji,jj) = 0.5_wp * vmask(ji,jj,1) * r1_e1e2v(ji,jj) & 781 & * ( e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji,jj+1) * sshn(ji,jj+1) ) 782 END DO 783 END DO 784 CALL lbc_lnk( zwv(:,:), 'V', 1._wp ) 785 z_tmax = MAXVAL( ssvmask(:,:) * ssvmask(:,:) * ABS( hv_0(:,:) + zwv(:,:) - zht(:,:) ) ) 520 786 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 521 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(e3t_b))) =', z_tmax 522 ! 523 z_tmax = MAXVAL( tmask(:,:,1) * ABS( sshb(:,:) ) ) 787 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(hv_0+sshv_n-SUM(e3v_n))) =', z_tmax 788 ! 789 zht(:,:) = 0.0_wp 790 DO jk = 1, jpkm1 791 DO jj = 1, jpjm1 792 zht(:,jj) = zht(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk)*umask(:,jj+1,jk) 793 END DO 794 END DO 795 zwu(:,:) = 0._wp 796 DO jj = 1, jpjm1 797 DO ji = 1, fs_jpim1 ! vector opt. 798 zwu(ji,jj) = 0.25_wp * umask(ji,jj,1) * umask(ji,jj+1,1) * r1_e1e2f(ji,jj) & 799 & * ( e1e2t(ji ,jj) * sshn(ji ,jj) + e1e2t(ji ,jj+1) * sshn(ji ,jj+1) & 800 & + e1e2t(ji+1,jj) * sshn(ji+1,jj) + e1e2t(ji+1,jj+1) * sshn(ji+1,jj+1) ) 801 END DO 802 END DO 803 CALL lbc_lnk( zht(:,:), 'F', 1._wp ) 804 CALL lbc_lnk( zwu(:,:), 'F', 1._wp ) 805 z_tmax = MAXVAL( fmask(:,:,1) * fmask(:,:,1) * ABS( hf_0(:,:) + zwu(:,:) - zht(:,:) ) ) 524 806 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 525 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(sshb))) =', z_tmax 526 ! 527 z_tmax = MAXVAL( tmask(:,:,1) * ABS( sshn(:,:) ) ) 528 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 529 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(sshn))) =', z_tmax 530 ! 531 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssha(:,:) ) ) 532 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 533 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssha))) =', z_tmax 807 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(hf_0+sshf_n-SUM(e3f_n))) =', z_tmax 808 ! 534 809 END IF 535 810 … … 699 974 !!---------------------------------------------------------------------- 700 975 ! 701 702 ! nmet = 1! Interface interpolation703 ! nmet = 2 ! Internal interfaces interpolation only, spread barotropic increment976 ! nmet = 0 ! Original method (Surely wrong) 977 nmet = 2 ! Interface interpolation 978 ! nmet = 2 ! Internal interfaces interpolation only, spread barotropic increment 704 979 ! Note that we kept surface weighted interpolation for barotropic increment to be compliant 705 980 ! with what is done in surface pressure module. … … 711 986 END IF 712 987 ! 713 ztap = 0. 0_wp ! Minimum fraction of T-point thickness at cell interfaces714 zsmall = 1.e- 3_wp ! Minimum thickness at U or V points (m)988 ztap = 0._wp ! Minimum fraction of T-point thickness at cell interfaces 989 zsmall = 1.e-8_wp ! Minimum thickness at U or V points (m) 715 990 ! 716 991 IF ( (nmet==1).OR.(nmet==2) ) THEN … … 744 1019 END DO 745 1020 ENDIF 1021 zw(:,:,:) = (zw(:,:,:) - gdepw_0(:,:,:))*tmask(:,:,:) 746 1022 ! 747 1023 END SELECT … … 768 1044 ! Correction at last level: 769 1045 jkbot = mbku(ji,jj) 770 zdo = hu_0(ji,jj)1046 zdo = 0._wp 771 1047 DO jk=jkbot,1,-1 772 zup = 0.5_wp * ( zw(ji ,jj,jk) + zw(ji+1,jj,jk) ) 1048 zup = 0.5_wp * ( e1e2t(ji ,jj)*zw(ji ,jj,jk) & 1049 & + e1e2t(ji+1,jj)*zw(ji+1,jj,jk) ) * r1_e1e2u(ji,jj) 773 1050 ! 774 1051 ! If there is a step, taper bottom interface: 775 IF ((hu_0(ji,jj) < 0.5_wp * ( ht_0(ji,jj) + ht_0(ji+1,jj) ) ).AND.(zup>zdo)) THEN776 IF ( ht_0(ji+1,jj) < ht_0(ji,jj) ) THEN777 zmin = ztap * (zw(ji+1,jj,jk+1)-zw(ji+1,jj,jk))778 ELSE779 zmin = ztap * (zw(ji ,jj,jk+1)-zw(ji ,jj,jk))780 ENDIF781 zup = MIN(zup, zdo-zmin)782 ENDIF783 zup = MIN(zup, zdo -zsmall)784 pe3_out(ji,jj,jk) = (zdo - zup - e3u_0(ji,jj,jk)) * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd )1052 ! IF ((hu_0(ji,jj) < 0.5_wp * ( ht_0(ji,jj) + ht_0(ji+1,jj) ) ).AND.(jk==jkbot)) THEN 1053 ! IF ( ht_0(ji+1,jj) < ht_0(ji,jj) ) THEN 1054 ! zmin = ztap * (zw(ji+1,jj,jk+1)-zw(ji+1,jj,jk)) 1055 ! ELSE 1056 ! zmin = ztap * (zw(ji ,jj,jk+1)-zw(ji ,jj,jk)) 1057 ! ENDIF 1058 ! zup = MIN(zup, zdo-zmin) 1059 ! ENDIF 1060 zup = MIN(zup, zdo+e3u_0(ji,jj,jk)-zsmall) 1061 pe3_out(ji,jj,jk) = (zdo - zup) * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) 785 1062 zdo = zup 786 1063 END DO … … 795 1072 pe3_out(ji,jj,jk) = pe3_out(ji,jj,jk) & 796 1073 & + ( pe3_out(ji,jj,jk) + e3u_0(ji,jj,jk) ) & 797 & / ( hu_0(ji,jj) + 1._wp - umask(ji,jj,1) ) &1074 & / ( hu_0(ji,jj) + 1._wp - ssumask(ji,jj) ) & 798 1075 & * 0.5_wp * r1_e1e2u(ji,jj) & 799 1076 & * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) & … … 824 1101 ! Correction at last level: 825 1102 jkbot = mbkv(ji,jj) 826 zdo = hv_0(ji,jj)1103 zdo = 0._wp 827 1104 DO jk=jkbot,1,-1 828 zup = 0.5_wp * ( zw(ji,jj ,jk) + zw(ji,jj+1,jk) ) 1105 zup = 0.5_wp * ( e1e2t(ji,jj ) * zw(ji,jj ,jk) & 1106 & + e1e2t(ji,jj+1) * zw(ji,jj+1,jk) ) * r1_e1e2v(ji,jj) 829 1107 ! 830 1108 ! If there is a step, taper bottom interface: 831 IF ((hv_0(ji,jj) < 0.5_wp * ( ht_0(ji,jj) + ht_0(ji,jj+1) ) ).AND.(zup>zdo)) THEN832 IF ( ht_0(ji,jj+1) < ht_0(ji,jj) ) THEN833 zmin = ztap * (zw(ji,jj+1,jk+1)-zw(ji,jj+1,jk))834 ELSE835 zmin = ztap * (zw(ji ,jj,jk+1)-zw(ji ,jj,jk))836 ENDIF837 zup = MIN(zup, zdo-zmin)838 ENDIF839 zup = MIN(zup, zdo -zsmall)840 pe3_out(ji,jj,jk) = (zdo - zup - e3v_0(ji,jj,jk)) * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd )1109 ! IF ((hv_0(ji,jj) < 0.5_wp * ( ht_0(ji,jj) + ht_0(ji,jj+1) ) ).AND.(jk==jkbot)) THEN 1110 ! IF ( ht_0(ji,jj+1) < ht_0(ji,jj) ) THEN 1111 ! zmin = ztap * (zw(ji,jj+1,jk+1)-zw(ji,jj+1,jk)) 1112 ! ELSE 1113 ! zmin = ztap * (zw(ji ,jj,jk+1)-zw(ji ,jj,jk)) 1114 ! ENDIF 1115 ! zup = MIN(zup, zdo-zmin) 1116 ! ENDIF 1117 zup = MIN(zup, zdo + e3v_0(ji,jj,jk) - zsmall) 1118 pe3_out(ji,jj,jk) = (zdo - zup) * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) 841 1119 zdo = zup 842 1120 END DO … … 851 1129 pe3_out(ji,jj,jk) = pe3_out(ji,jj,jk) & 852 1130 & + ( pe3_out(ji,jj,jk) + e3v_0(ji,jj,jk) ) & 853 & / ( hv_0(ji,jj) + 1._wp - vmask(ji,jj,1) ) &1131 & / ( hv_0(ji,jj) + 1._wp - ssvmask(ji,jj) ) & 854 1132 & * 0.5_wp * r1_e1e2v(ji,jj) & 855 1133 * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) & … … 883 1161 ! bottom correction: 884 1162 jkbot = MIN(mbku(ji,jj), mbku(ji,jj+1)) 885 zdo = hf_0(ji,jj)1163 zdo = 0._wp 886 1164 DO jk=jkbot,1,-1 887 zup = 0.25_wp * ( zw(ji ,jj ,jk) &888 & + zw(ji+1,jj ,jk) &889 & + zw(ji ,jj+1,jk) &890 & + zw(ji+1,jj+1,jk))1165 zup = 0.25_wp * ( e1e2t(ji ,jj ) * zw(ji ,jj ,jk) & 1166 & + e1e2t(ji+1,jj ) * zw(ji+1,jj ,jk) & 1167 & + e1e2t(ji ,jj+1) * zw(ji ,jj+1,jk) & 1168 & + e1e2t(ji+1,jj+1) * zw(ji+1,jj+1,jk) ) * r1_e1e2f(ji,jj) 891 1169 ! 892 1170 ! If there is a step, taper bottom interface: 893 IF ((hf_0(ji,jj) < 0.5_wp * ( hu_0(ji,jj ) + hu_0(ji,jj+1) ) ).AND.(zup>zdo)) THEN894 IF ( hu_0(ji,jj+1) < hu_0(ji,jj) ) THEN895 IF ( ht_0(ji+1,jj+1) < ht_0(ji ,jj+1) ) THEN896 zmin = ztap * (zw(ji+1,jj+1,jk+1)-zw(ji+1,jj+1,jk))897 ELSE898 zmin = ztap * (zw(ji ,jj+1,jk+1)-zw(ji ,jj+1,jk))899 ENDIF900 ELSE901 IF ( ht_0(ji+1,jj ) < ht_0(ji ,jj ) ) THEN902 zmin = ztap * (zw(ji+1,jj ,jk+1)-zw(ji+1,jj ,jk))903 ELSE904 zmin = ztap * (zw(ji ,jj ,jk+1)-zw(ji ,jj ,jk))905 ENDIF906 ENDIF907 zup = MIN(zup, zdo-zmin)908 ENDIF909 zup = MIN(zup, zdo -zsmall)1171 ! IF ((hf_0(ji,jj) < 0.5_wp * ( hu_0(ji,jj ) + hu_0(ji,jj+1) ) ).AND.(jk==jkbot)) THEN 1172 ! IF ( hu_0(ji,jj+1) < hu_0(ji,jj) ) THEN 1173 ! IF ( ht_0(ji+1,jj+1) < ht_0(ji ,jj+1) ) THEN 1174 ! zmin = ztap * (zw(ji+1,jj+1,jk+1)-zw(ji+1,jj+1,jk)) 1175 ! ELSE 1176 ! zmin = ztap * (zw(ji ,jj+1,jk+1)-zw(ji ,jj+1,jk)) 1177 ! ENDIF 1178 ! ELSE 1179 ! IF ( ht_0(ji+1,jj ) < ht_0(ji ,jj ) ) THEN 1180 ! zmin = ztap * (zw(ji+1,jj ,jk+1)-zw(ji+1,jj ,jk)) 1181 ! ELSE 1182 ! zmin = ztap * (zw(ji ,jj ,jk+1)-zw(ji ,jj ,jk)) 1183 ! ENDIF 1184 ! ENDIF 1185 ! zup = MIN(zup, zdo-zmin) 1186 ! ENDIF 1187 zup = MIN(zup, zdo+e3f_0(ji,jj,jk)-zsmall) 910 1188 ! 911 pe3_out(ji,jj,jk) = ( zdo - zup - e3f_0(ji,jj,jk)) &1189 pe3_out(ji,jj,jk) = ( zdo - zup ) & 912 1190 & *( umask(ji,jj,jk) * umask(ji,jj+1,jk) * (1.0_wp - zlnwd) + zlnwd ) 913 1191 zdo = zup … … 1006 1284 id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. ) 1007 1285 id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) 1008 id5 = iom_varid( numror, 'hdiv _lf', ldstop = .FALSE. )1286 id5 = iom_varid( numror, 'hdivn_lf', ldstop = .FALSE. ) 1009 1287 ! ! --------- ! 1010 1288 ! ! all cases ! … … 1068 1346 ! ! ------------ ! 1069 1347 IF( id5 > 0 ) THEN ! required array exists 1070 CALL iom_get( numror, jpdom_autoglo, 'hdiv _lf', hdiv_lf(:,:,:), ldxios = lrxios )1348 CALL iom_get( numror, jpdom_autoglo, 'hdivn_lf', hdivn_lf(:,:,:), ldxios = lrxios ) 1071 1349 ELSE ! array is missing 1072 hdiv _lf(:,:,:) = 0.0_wp1350 hdivn_lf(:,:,:) = 0.0_wp 1073 1351 ENDIF 1074 1352 ENDIF … … 1140 1418 tilde_e3t_b(:,:,:) = 0._wp 1141 1419 tilde_e3t_n(:,:,:) = 0._wp 1142 IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0._wp 1420 IF( ln_vvl_ztilde ) THEN 1421 hdivn_lf(:,:,:) = 0._wp 1422 un_lf(:,:,:) = 0._wp 1423 vn_lf(:,:,:) = 0._wp 1424 ENDIF 1143 1425 END IF 1144 1426 ENDIF … … 1162 1444 IF( ln_vvl_ztilde ) THEN ! z_tilde case ! 1163 1445 ! ! ------------ ! 1164 CALL iom_rstput( kt, nitrst, numrow, 'hdiv _lf', hdiv_lf(:,:,:), ldxios = lwxios)1446 CALL iom_rstput( kt, nitrst, numrow, 'hdivn_lf', hdivn_lf(:,:,:), ldxios = lwxios) 1165 1447 ENDIF 1166 1448 ! … … 1179 1461 !!---------------------------------------------------------------------- 1180 1462 INTEGER :: ioptio, ios 1181 !! 1182 NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, & 1183 & ln_vvl_zstar_at_eqtor , rn_ahe3 , rn_rst_e3t , & 1184 & rn_lf_cutoff , rn_zdef_max , ln_vvl_dbg ! not yet implemented: ln_vvl_kepe 1185 !!---------------------------------------------------------------------- 1463 1464 NAMELIST/nam_vvl/ ln_vvl_zstar , ln_vvl_ztilde , & 1465 & ln_vvl_layer , ln_vvl_ztilde_as_zstar , & 1466 & ln_vvl_zstar_at_eqtor , ln_vvl_zstar_on_shelf , & 1467 & ln_vvl_adv_cn2 , ln_vvl_adv_fct , & 1468 & ln_vvl_lap , ln_vvl_blp , & 1469 & rn_ahe3_lap , rn_ahe3_blp , & 1470 & rn_rst_e3t , rn_lf_cutoff , & 1471 & ln_vvl_regrid , & 1472 & ln_vvl_ramp , rn_day_ramp , & 1473 & ln_vvl_dbg ! not yet implemented: ln_vvl_kepe 1474 !!---------------------------------------------------------------------- 1186 1475 ! 1187 1476 REWIND( numnam_ref ) ! Namelist nam_vvl in reference namelist : … … 1197 1486 WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate' 1198 1487 WRITE(numout,*) '~~~~~~~~~~~' 1199 WRITE(numout,*) ' Namelist nam_vvl : chose a vertical coordinate' 1200 WRITE(numout,*) ' zstar ln_vvl_zstar = ', ln_vvl_zstar 1201 WRITE(numout,*) ' ztilde ln_vvl_ztilde = ', ln_vvl_ztilde 1202 WRITE(numout,*) ' layer ln_vvl_layer = ', ln_vvl_layer 1203 WRITE(numout,*) ' ztilde as zstar ln_vvl_ztilde_as_zstar = ', ln_vvl_ztilde_as_zstar 1488 WRITE(numout,*) ' Namelist nam_vvl : chose a vertical coordinate' 1489 WRITE(numout,*) ' zstar ln_vvl_zstar = ', ln_vvl_zstar 1490 WRITE(numout,*) ' ztilde ln_vvl_ztilde = ', ln_vvl_ztilde 1491 WRITE(numout,*) ' layer ln_vvl_layer = ', ln_vvl_layer 1492 WRITE(numout,*) ' ztilde as zstar ln_vvl_ztilde_as_zstar = ', ln_vvl_ztilde_as_zstar 1493 ! WRITE(numout,*) ' Namelist nam_vvl : chose kinetic-to-potential energy conservation' 1494 ! WRITE(numout,*) ' ln_vvl_kepe = ', ln_vvl_kepe 1204 1495 WRITE(numout,*) ' ztilde near the equator ln_vvl_zstar_at_eqtor = ', ln_vvl_zstar_at_eqtor 1205 WRITE(numout,*) ' !' 1206 WRITE(numout,*) ' thickness diffusion coefficient rn_ahe3 = ', rn_ahe3 1207 WRITE(numout,*) ' maximum e3t deformation fractional change rn_zdef_max = ', rn_zdef_max 1496 WRITE(numout,*) ' ztilde on shelves ln_vvl_zstar_on_shelf = ', ln_vvl_zstar_on_shelf 1497 WRITE(numout,*) ' Namelist nam_vvl : thickness advection scheme' 1498 WRITE(numout,*) ' 2nd order ln_vvl_adv_cn2 = ', ln_vvl_adv_cn2 1499 WRITE(numout,*) ' 2nd order FCT ln_vvl_adv_fct = ', ln_vvl_adv_fct 1500 WRITE(numout,*) ' Namelist nam_vvl : thickness diffusion scheme' 1501 WRITE(numout,*) ' Laplacian ln_vvl_lap = ', ln_vvl_lap 1502 WRITE(numout,*) ' Bilaplacian ln_vvl_blp = ', ln_vvl_blp 1503 WRITE(numout,*) ' Laplacian coefficient rn_ahe3_lap = ', rn_ahe3_lap 1504 WRITE(numout,*) ' Bilaplacian coefficient rn_ahe3_blp = ', rn_ahe3_blp 1505 WRITE(numout,*) ' Namelist nam_vvl : layers regriding' 1506 WRITE(numout,*) ' ln_vvl_regrid = ', ln_vvl_regrid 1507 WRITE(numout,*) ' Namelist nam_vvl : linear ramp at startup' 1508 WRITE(numout,*) ' ln_vvl_ramp = ', ln_vvl_ramp 1509 WRITE(numout,*) ' rn_day_ramp = ', rn_day_ramp 1208 1510 IF( ln_vvl_ztilde_as_zstar ) THEN 1209 WRITE(numout,*) ' ztilde running in zstar emulation mode (ln_vvl_ztilde_as_zstar=T)'1210 WRITE(numout,*) ' ignoring namelist timescale parameters and using:'1211 WRITE(numout,*) ' hard-wired : z-tilde to zstar restoration timescale (days)'1212 WRITE(numout,*) ' rn_rst_e3t = 0.e0'1213 WRITE(numout,*) ' hard-wired : z-tilde cutoff frequency of low-pass filter (days)'1214 WRITE(numout,*) ' rn_lf_cutoff =1.0/rdt'1511 WRITE(numout,*) ' ztilde running in zstar emulation mode; ' 1512 WRITE(numout,*) ' ignoring namelist timescale parameters and using:' 1513 WRITE(numout,*) ' hard-wired : z-tilde to zstar restoration timescale (days)' 1514 WRITE(numout,*) ' rn_rst_e3t = 0.0' 1515 WRITE(numout,*) ' hard-wired : z-tilde cutoff frequency of low-pass filter (days)' 1516 WRITE(numout,*) ' rn_lf_cutoff = 1.0/rdt' 1215 1517 ELSE 1216 WRITE(numout,*) ' z-tilde to zstar restoration timescale (days) rn_rst_e3t = ', rn_rst_e3t 1217 WRITE(numout,*) ' z-tilde cutoff frequency of low-pass filter (days) rn_lf_cutoff = ', rn_lf_cutoff 1518 WRITE(numout,*) ' Namelist nam_vvl : z-tilde to zstar restoration timescale (days)' 1519 WRITE(numout,*) ' rn_rst_e3t = ', rn_rst_e3t 1520 WRITE(numout,*) ' Namelist nam_vvl : z-tilde cutoff frequency of low-pass filter (days)' 1521 WRITE(numout,*) ' rn_lf_cutoff = ', rn_lf_cutoff 1218 1522 ENDIF 1219 WRITE(numout,*) ' debug prints flag ln_vvl_dbg = ', ln_vvl_dbg 1523 WRITE(numout,*) ' Namelist nam_vvl : debug prints' 1524 WRITE(numout,*) ' ln_vvl_dbg = ', ln_vvl_dbg 1525 ENDIF 1526 ! 1527 IF ( ln_vvl_ztilde.OR.ln_vvl_layer ) THEN 1528 ioptio = 0 ! Choose one advection scheme at most 1529 IF( ln_vvl_adv_cn2 ) ioptio = ioptio + 1 1530 IF( ln_vvl_adv_fct ) ioptio = ioptio + 1 1531 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE thickness advection scheme in namelist nam_vvl' ) 1220 1532 ENDIF 1221 1533 ! … … 1237 1549 ENDIF 1238 1550 ! 1551 ! Use of "shelf horizon depths" should be allowed with s-z coordinates, but we restrict it to zco and zps 1552 ! for the time being 1553 IF ( ln_sco ) THEN 1554 ll_shorizd=.FALSE. 1555 ELSE 1556 ll_shorizd=.TRUE. 1557 ENDIF 1558 ! 1239 1559 #if defined key_agrif 1240 1560 IF( (.NOT.Agrif_Root()).AND.(.NOT.ln_vvl_zstar) ) CALL ctl_stop( 'AGRIF is implemented with zstar coordinate only' ) … … 1243 1563 END SUBROUTINE dom_vvl_ctl 1244 1564 1565 SUBROUTINE dom_vvl_regrid( kt ) 1566 !!---------------------------------------------------------------------- 1567 !! *** ROUTINE dom_vvl_regrid *** 1568 !! 1569 !! ** Purpose : Ensure "well-behaved" vertical grid 1570 !! 1571 !! ** Method : More or less adapted from references below. 1572 !!regrid 1573 !! ** Action : Ensure that thickness are above a given value, spaced enough 1574 !! and revert to Eulerian coordinates near the bottom. 1575 !! 1576 !! References : Bleck, R. and S. Benjamin, 1993: Regional Weather Prediction 1577 !! with a Model Combining Terrain-following and Isentropic 1578 !! coordinates. Part I: Model Description. Monthly Weather Rev., 1579 !! 121, 1770-1785. 1580 !! Toy, M., 2011: Incorporating Condensational Heating into a 1581 !! Nonhydrostatic Atmospheric Model Based on a Hybrid Isentropic- 1582 !! Sigma Vertical Coordinate. Monthly Weather Rev., 139, 2940-2954. 1583 !!---------------------------------------------------------------------- 1584 !! * Arguments 1585 INTEGER, INTENT( in ) :: kt ! time step 1586 1587 !! * Local declarations 1588 INTEGER :: ji, jj, jk ! dummy loop indices 1589 LOGICAL :: ll_chk_bot2top, ll_chk_top2bot, ll_lapdiff_cond 1590 LOGICAL :: ll_zdiff_cond, ll_blpdiff_cond 1591 INTEGER :: jkbot 1592 REAL(wp) :: zh_min, zh_0, zh2, zdiff, zh_max, ztmph, ztmpd 1593 REAL(wp) :: zufim1, zufi, zvfjm1, zvfj, dzmin_int, dzmin_surf 1594 REAL(wp) :: zh_new, zh_old, zh_bef, ztmp, ztmp1, z2dt, zh_up, zh_dwn 1595 REAL(wp) :: zeu2, zev2, zfrch_stp, zfrch_rel, zfrac_bot, zscal_bot 1596 REAL(wp) :: zhdiff, zhdiff2, zvdiff, zhlim, zhlim2, zvlim 1597 REAL(wp), DIMENSION(jpi,jpj) :: zdw, zwu, zwv 1598 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwdw, zwdw_b 1599 !!---------------------------------------------------------------------- 1600 1601 IF( ln_timing ) CALL timing_start('dom_vvl_regrid') 1602 ! 1603 ! 1604 ! Some user defined parameters below: 1605 ll_chk_bot2top = .TRUE. 1606 ll_chk_top2bot = .TRUE. 1607 dzmin_int = 0.1_wp ! Absolute minimum depth in the interior (in meters) 1608 dzmin_surf = 1.0_wp ! Absolute minimum depth at the surface (in meters) 1609 zfrch_stp = 0.5_wp ! Maximum fractionnal thickness change in one time step (<= 1.) 1610 zfrch_rel = 0.4_wp ! Maximum relative thickness change in the vertical (<= 1.) 1611 zfrac_bot = 0.05_wp ! Fraction of bottom level allowed to change 1612 zscal_bot = 2.0_wp ! Depth lengthscale 1613 ll_zdiff_cond = .TRUE. ! Conditionnal vertical diffusion of interfaces 1614 zvdiff = 0.2_wp ! m 1615 zvlim = 0.5_wp ! max d2h/dh 1616 ll_lapdiff_cond = .TRUE. ! Conditionnal Laplacian diffusion of interfaces 1617 zhdiff = 0.01_wp ! ad. 1618 zhlim = 0.03_wp ! ad. max lap(z)*e1 1619 ll_blpdiff_cond = .FALSE. ! Conditionnal Bilaplacian diffusion of interfaces 1620 zhdiff2 = 0.2_wp ! ad. 1621 ! zhlim2 = 3.e-11_wp ! max bilap(z) 1622 zhlim2 = 0.01_wp ! ad. max bilap(z)*e1**3 1623 ! --------------------------------------------------------------------------------------- 1624 ! 1625 ! Set arrays determining maximum vertical displacement at the bottom: 1626 !-------------------------------------------------------------------- 1627 IF ( kt==nit000 ) THEN 1628 DO jj = 2, jpjm1 1629 DO ji = 2, jpim1 1630 jk = MIN(mbkt(ji,jj), mbkt(ji+1,jj), mbkt(ji-1,jj), mbkt(ji,jj+1), mbkt(ji,jj-1)) 1631 jk = MIN(jk,mbkt(ji-1,jj-1), mbkt(ji-1,jj+1), mbkt(ji+1,jj+1), mbkt(ji+1,jj-1)) 1632 i_int_bot(ji,jj) = jk 1633 END DO 1634 END DO 1635 dsm(:,:) = REAL( i_int_bot(:,:), wp ) ; CALL lbc_lnk(dsm(:,:),'T',1.) 1636 i_int_bot(:,:) = MAX( INT( dsm(:,:) ), 1 ) 1637 1638 DO jj = 2, jpjm1 1639 DO ji = 2, jpim1 1640 zdw(ji,jj) = MAX(ABS(ht_0(ji,jj)-ht_0(ji+1,jj))*umask(ji ,jj,1), & 1641 & ABS(ht_0(ji,jj)-ht_0(ji-1,jj))*umask(ji-1,jj,1), & 1642 & ABS(ht_0(ji,jj)-ht_0(ji,jj+1))*vmask(ji,jj ,1), & 1643 & ABS(ht_0(ji,jj)-ht_0(ji,jj-1))*vmask(ji,jj-1,1) ) 1644 zdw(ji,jj) = MAX(zscal_bot * zdw(ji,jj), rsmall ) 1645 END DO 1646 END DO 1647 CALL lbc_lnk( zdw(:,:), 'T', 1. ) 1648 1649 DO jj = 2, jpjm1 1650 DO ji = 2, jpim1 1651 dsm(ji,jj) = 1._wp/16._wp * ( zdw(ji-1,jj-1) + zdw(ji+1,jj-1) & 1652 & + zdw(ji-1,jj+1) + zdw(ji+1,jj+1) & 1653 & + 2._wp*( zdw(ji ,jj-1) + zdw(ji-1,jj ) & 1654 & + zdw(ji+1,jj ) + zdw(ji ,jj+1) ) & 1655 & + 4._wp* zdw(ji ,jj ) ) 1656 END DO 1657 END DO 1658 1659 CALL lbc_lnk( dsm(:,:), 'T', 1. ) 1660 1661 IF (ln_zps) THEN 1662 DO jj = 1, jpj 1663 DO ji = 1, jpi 1664 jk = i_int_bot(ji,jj) 1665 hsm(ji,jj) = zfrac_bot * e3w_1d(jk) 1666 ! dsm(ji,jj) = MAX(dsm(ji,jj), 0.1_wp*ht_0(ji,jj)) 1667 END DO 1668 END DO 1669 ELSE 1670 DO jj = 1, jpj 1671 DO ji = 1, jpi 1672 jk = i_int_bot(ji,jj) 1673 hsm(ji,jj) = zfrac_bot * e3w_0(ji,jj,jk) 1674 ! dsm(ji,jj) = MAX(dsm(ji,jj), 0.1_wp*ht_0(ji,jj)) 1675 END DO 1676 END DO 1677 ENDIF 1678 END IF 1679 1680 ! Provisionnal interface depths: 1681 !------------------------------- 1682 zwdw(:,:,1) = 0.e0 1683 DO jj = 1, jpj 1684 DO ji = 1, jpi 1685 DO jk = 2, jpk 1686 zwdw(ji,jj,jk) = zwdw(ji,jj,jk-1) + & 1687 & (tilde_e3t_a(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1) 1688 END DO 1689 END DO 1690 END DO 1691 ! 1692 ! Conditionnal horizontal Laplacian diffusion: 1693 !--------------------------------------------- 1694 IF ( ll_lapdiff_cond ) THEN 1695 ! 1696 zwdw_b(:,:,1) = 0._wp 1697 DO jj = 1, jpj 1698 DO ji = 1, jpi 1699 DO jk=2,jpk 1700 zwdw_b(ji,jj,jk) = zwdw_b(ji,jj,jk-1) + & 1701 & (tilde_e3t_b(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1) 1702 END DO 1703 END DO 1704 END DO 1705 ! 1706 DO jk = 2, jpkm1 1707 zwu(:,:) = 0._wp 1708 zwv(:,:) = 0._wp 1709 DO jj = 1, jpjm1 1710 DO ji = 1, fs_jpim1 ! vector opt. 1711 zwu(ji,jj) = umask(ji,jj,jk) * e2_e1u(ji,jj) & 1712 & * ( zwdw_b(ji,jj,jk) - zwdw_b(ji+1,jj ,jk) ) 1713 zwv(ji,jj) = vmask(ji,jj,jk) * e1_e2v(ji,jj) & 1714 & * ( zwdw_b(ji,jj,jk) - zwdw_b(ji ,jj+1,jk) ) 1715 END DO 1716 END DO 1717 DO jj = 2, jpjm1 1718 DO ji = fs_2, fs_jpim1 ! vector opt. 1719 ztmp1 = ( zwu(ji-1,jj ) - zwu(ji,jj) & 1720 & + zwv(ji ,jj-1) - zwv(ji,jj) ) * r1_e1e2t(ji,jj) 1721 zh2 = MAX(abs(ztmp1)-zhlim*SQRT(r1_e1e2t(ji,jj)), 0._wp) 1722 ztmp = SIGN(zh2, ztmp1) 1723 zeu2 = zhdiff * e1e2t(ji,jj)*e1e2t(ji,jj)/(e1t(ji,jj)*e1t(ji,jj) + e2t(ji,jj)*e2t(ji,jj)) 1724 zwdw(ji,jj,jk) = zwdw(ji,jj,jk) + zeu2 * ztmp * tmask(ji,jj,jk) 1725 END DO 1726 END DO 1727 END DO 1728 ! 1729 ENDIF 1730 1731 ! Conditionnal horizontal Bilaplacian diffusion: 1732 !----------------------------------------------- 1733 IF ( ll_blpdiff_cond ) THEN 1734 ! 1735 zwdw_b(:,:,1) = 0._wp 1736 DO jj = 1, jpj 1737 DO ji = 1, jpi 1738 DO jk = 2,jpkm1 1739 zwdw_b(ji,jj,jk) = zwdw_b(ji,jj,jk-1) + & 1740 & (tilde_e3t_b(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1) 1741 END DO 1742 END DO 1743 END DO 1744 ! 1745 DO jk = 2, jpkm1 1746 zwu(:,:) = 0._wp 1747 zwv(:,:) = 0._wp 1748 DO jj = 1, jpjm1 1749 DO ji = 1, fs_jpim1 ! vector opt. 1750 zwu(ji,jj) = umask(ji,jj,jk) * e2_e1u(ji,jj) & 1751 & * ( zwdw_b(ji,jj,jk) - zwdw_b(ji+1,jj ,jk) ) 1752 zwv(ji,jj) = vmask(ji,jj,jk) * e1_e2v(ji,jj) & 1753 & * ( zwdw_b(ji,jj,jk) - zwdw_b(ji ,jj+1,jk) ) 1754 END DO 1755 END DO 1756 DO jj = 2, jpjm1 1757 DO ji = fs_2, fs_jpim1 ! vector opt. 1758 zwdw_b(ji,jj,jk) = -( (zwu(ji-1,jj ) - zwu(ji,jj)) & 1759 & + (zwv(ji ,jj-1) - zwv(ji,jj)) ) * r1_e1e2t(ji,jj) 1760 END DO 1761 END DO 1762 END DO 1763 ! 1764 CALL lbc_lnk( zwdw_b(:,:,:), 'T', 1. ) 1765 ! 1766 DO jk = 2, jpkm1 1767 zwu(:,:) = 0._wp 1768 zwv(:,:) = 0._wp 1769 DO jj = 1, jpjm1 1770 DO ji = 1, fs_jpim1 ! vector opt. 1771 zwu(ji,jj) = umask(ji,jj,jk) * e2_e1u(ji,jj) & 1772 & * ( zwdw_b(ji,jj,jk) - zwdw_b(ji+1,jj ,jk) ) 1773 zwv(ji,jj) = vmask(ji,jj,jk) * e1_e2v(ji,jj) & 1774 & * ( zwdw_b(ji,jj,jk) - zwdw_b(ji ,jj+1,jk) ) 1775 END DO 1776 END DO 1777 DO jj = 2, jpjm1 1778 DO ji = fs_2, fs_jpim1 ! vector opt. 1779 ztmp1 = ( (zwu(ji-1,jj ) - zwv(ji,jj)) & 1780 & + (zwv(ji ,jj-1) - zwv(ji,jj)) ) * r1_e1e2t(ji,jj) 1781 zh2 = MAX(abs(ztmp1)-zhlim2*SQRT(r1_e1e2t(ji,jj))*r1_e1e2t(ji,jj), 0._wp) 1782 ztmp = SIGN(zh2, ztmp1) 1783 zeu2 = zhdiff2 * e1e2t(ji,jj)*e1e2t(ji,jj) / 16._wp 1784 zwdw(ji,jj,jk) = zwdw(ji,jj,jk) + zeu2 * ztmp * tmask(ji,jj,jk) 1785 END DO 1786 END DO 1787 END DO 1788 ! 1789 ENDIF 1790 1791 ! Conditionnal vertical diffusion: 1792 !--------------------------------- 1793 IF ( ll_zdiff_cond ) THEN 1794 DO jk = 2, jpkm1 1795 DO jj = 2, jpjm1 1796 DO ji = fs_2, fs_jpim1 ! vector opt. 1797 ztmp = -( (tilde_e3t_b(ji,jj,jk-1)+e3t_0(ji,jj,jk-1))*tmask(ji,jj,jk-1) & 1798 -(tilde_e3t_b(ji,jj,jk )+e3t_0(ji,jj,jk ))*tmask(ji,jj,jk ) ) 1799 ztmp1 = 0.5_wp * ( tilde_e3t_b(ji,jj,jk-1) + e3t_0(ji,jj,jk-1) & 1800 & +tilde_e3t_b(ji,jj,jk ) + e3t_0(ji,jj,jk ) ) 1801 zh2 = MAX(abs(ztmp)-zvlim*ztmp1, 0._wp) 1802 ztmp = SIGN(zh2, ztmp) 1803 IF ((jk==mbkt(ji,jj)).AND.(ln_zps)) ztmp=0.e0 1804 zwdw(ji,jj,jk) = zwdw(ji,jj,jk) + zvdiff * ztmp * tmask(ji,jj,jk) 1805 END DO 1806 END DO 1807 END DO 1808 ENDIF 1809 ! 1810 ! Check grid from the bottom to the surface 1811 !------------------------------------------ 1812 IF ( ll_chk_bot2top ) THEN 1813 DO jj = 2, jpjm1 1814 DO ji = 2, jpim1 1815 jkbot = mbkt(ji,jj) 1816 DO jk = jkbot,2,-1 1817 ! 1818 zh_0 = e3t_0(ji,jj,jk) 1819 zh_bef = MIN(tilde_e3t_b(ji,jj,jk) + zh_0, tilde_e3t_b(ji,jj,jk-1) + e3t_0(ji,jj,jk-1)) 1820 zh_old = zwdw(ji,jj,jk+1) - zwdw(ji,jj,jk) 1821 zh_min = MIN(zh_0/3._wp, dzmin_int) 1822 ! 1823 ! Set maximum and minimum vertical excursions 1824 ztmph = hsm(ji,jj) 1825 ztmpd = dsm(ji,jj) 1826 zh2 = ztmph * exp(-(gdepw_0(ji,jj,jk)-gdepw_0(ji,jj,i_int_bot(ji,jj)))/ztmpd) 1827 ! zh2 = ztmph * exp(-(gdepw_0(ji,jj,jk)-gdepw_0(ji,jj,i_int_bot(ji,jj)+1))/ztmpd) 1828 zh2 = MAX(zh2,0.001_wp) ! Extend tolerance a bit for stability reasons (to be explored) 1829 zdiff = cush_max(gdepw_0(ji,jj,jk)-zwdw(ji,jj,jk), zh2 ) 1830 zwdw(ji,jj,jk) = MAX(zwdw(ji,jj,jk), gdepw_0(ji,jj,jk) - zdiff) 1831 zdiff = cush_max(zwdw(ji,jj,jk)-gdepw_0(ji,jj,jk), zh2 ) 1832 zwdw(ji,jj,jk) = MIN(zwdw(ji,jj,jk), gdepw_0(ji,jj,jk) + zdiff) 1833 ! 1834 ! New layer thickness: 1835 zh_new = zwdw(ji,jj,jk+1) - zwdw(ji,jj,jk) 1836 ! 1837 ! Ensure minimum layer thickness: 1838 ! zh_new = MAX((1._wp-zfrch_stp)*zh_bef, zh_new) 1839 zh_new = cush(zh_new, zh_min) 1840 ! 1841 ! Final flux: 1842 zdiff = (zh_new - zh_old) * tmask(ji,jj,jk) 1843 ! 1844 ! Limit thickness change in 1 time step: 1845 ztmp = MIN( ABS(zdiff), zfrch_stp*zh_bef ) 1846 zdiff = SIGN(ztmp, zh_new - zh_old) 1847 zh_new = zdiff + zh_old 1848 ! 1849 zwdw(ji,jj,jk) = zwdw(ji,jj,jk+1) - zh_new 1850 END DO 1851 END DO 1852 END DO 1853 END IF 1854 ! 1855 ! Check grid from the surface to the bottom 1856 !------------------------------------------ 1857 IF ( ll_chk_top2bot ) THEN 1858 DO jj = 2, jpjm1 1859 DO ji = 2, jpim1 1860 jkbot = mbkt(ji,jj) 1861 DO jk = 1, jkbot-1 1862 ! 1863 zh_0 = e3t_0(ji,jj,jk) 1864 zh_bef = MIN(tilde_e3t_b(ji,jj,jk) + zh_0, tilde_e3t_b(ji,jj,jk+1) + e3t_0(ji,jj,jk+1)) 1865 zh_old = zwdw(ji,jj,jk+1) - zwdw(ji,jj,jk) 1866 zh_min = MIN(zh_0/3._wp, dzmin_int) 1867 ! 1868 zwdw(ji,jj,jk+1) = MAX(zwdw(ji,jj,jk+1), REAL(jk)*dzmin_surf) 1869 ! 1870 ! New layer thickness: 1871 zh_new = zwdw(ji,jj,jk+1) - zwdw(ji,jj,jk) 1872 ! 1873 ! Ensure minimum layer thickness: 1874 ! zh_new = MAX((1._wp-zfrch_stp)*zh_bef, zh_new) 1875 zh_new = cush(zh_new, zh_min) 1876 ! 1877 ! Final flux: 1878 zdiff = (zh_new -zh_old) * tmask(ji,jj,jk) 1879 ! 1880 ! Limit flux: 1881 ! ztmp = MIN( ABS(zdiff), zfrch_stp*zh_bef ) 1882 ! zdiff = SIGN(ztmp, zh_new - zh_old) 1883 zh_new = zdiff + zh_old 1884 ! 1885 zwdw(ji,jj,jk+1) = zwdw(ji,jj,jk) + zh_new 1886 END DO 1887 ! 1888 END DO 1889 END DO 1890 ENDIF 1891 ! 1892 DO jj = 2, jpjm1 1893 DO ji = 2, jpim1 1894 DO jk = 1, jpkm1 1895 tilde_e3t_a(ji,jj,jk) = (zwdw(ji,jj,jk+1)-zwdw(ji,jj,jk)-e3t_0(ji,jj,jk)) * tmask(ji,jj,jk) 1896 END DO 1897 END DO 1898 END DO 1899 ! 1900 ! 1901 IF( ln_timing ) CALL timing_stop('dom_vvl_regrid') 1902 ! 1903 END SUBROUTINE dom_vvl_regrid 1904 1905 FUNCTION cush(hin, hmin) RESULT(hout) 1906 !!---------------------------------------------------------------------- 1907 !! *** FUNCTION cush *** 1908 !! 1909 !! ** Purpose : 1910 !! 1911 !! ** Method : 1912 !! 1913 !!---------------------------------------------------------------------- 1914 IMPLICIT NONE 1915 REAL(wp), INTENT(in) :: hin, hmin 1916 REAL(wp) :: hout, zx, zh_cri 1917 !!---------------------------------------------------------------------- 1918 zh_cri = 3._wp * hmin 1919 ! 1920 IF ( hin<=0._wp ) THEN 1921 hout = hmin 1922 ! 1923 ELSEIF ( (hin>0._wp).AND.(hin<=zh_cri) ) THEN 1924 zx = hin/zh_cri 1925 hout = hmin * (1._wp + zx + zx*zx) 1926 ! 1927 ELSEIF ( hin>zh_cri ) THEN 1928 hout = hin 1929 ! 1930 ENDIF 1931 ! 1932 END FUNCTION cush 1933 1934 FUNCTION cush_max(hin, hmax) RESULT(hout) 1935 !!---------------------------------------------------------------------- 1936 !! *** FUNCTION cush *** 1937 !! 1938 !! ** Purpose : 1939 !! 1940 !! ** Method : 1941 !! 1942 !!---------------------------------------------------------------------- 1943 IMPLICIT NONE 1944 REAL(wp), INTENT(in) :: hin, hmax 1945 REAL(wp) :: hout, hmin, zx, zh_cri 1946 !!---------------------------------------------------------------------- 1947 hmin = 0.1_wp * hmax 1948 zh_cri = 3._wp * hmin 1949 ! 1950 IF ( (hin>=(hmax-zh_cri)).AND.(hin<=(hmax-hmin))) THEN 1951 zx = (hmax-hin)/zh_cri 1952 hout = hmax - hmin * (1._wp + zx + zx*zx) 1953 ! 1954 ELSEIF ( hin>(hmax-zh_cri) ) THEN 1955 hout = hmax - hmin 1956 ! 1957 ELSE 1958 hout = hin 1959 ! 1960 ENDIF 1961 ! 1962 END FUNCTION cush_max 1963 1964 SUBROUTINE dom_vvl_adv_fct( kt, pta, uin, vin ) 1965 !!---------------------------------------------------------------------- 1966 !! *** ROUTINE dom_vvl_adv_fct *** 1967 !! 1968 !! ** Purpose : Do thickness advection 1969 !! 1970 !! ** Method : FCT scheme to ensure positivity 1971 !! 1972 !! ** Action : - Update pta thickness tendency and diffusive fluxes 1973 !! - this is the total trend, hence it does include sea level motions 1974 !!---------------------------------------------------------------------- 1975 ! 1976 INTEGER, INTENT( in ) :: kt ! ocean time-step index 1977 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pta ! thickness baroclinic trend 1978 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: uin, vin ! input velocities 1979 ! 1980 INTEGER :: ji, jj, jk, ib, ib_bdy ! dummy loop indices 1981 INTEGER :: ikbu, ikbv, ibot 1982 REAL(wp) :: z2dtt, zbtr, ztra ! local scalar 1983 REAL(wp) :: zdi, zdj, zmin ! - - 1984 REAL(wp) :: zfp_ui, zfp_vj ! - - 1985 REAL(wp) :: zfm_ui, zfm_vj ! - - 1986 REAL(wp) :: zfp_hi, zfp_hj ! - - 1987 REAL(wp) :: zfm_hi, zfm_hj ! - - 1988 REAL(wp) :: ztout, ztin, zfac ! - - 1989 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwx, zwy, zwi 1990 !!---------------------------------------------------------------------- 1991 ! 1992 IF( ln_timing ) CALL timing_start('dom_vvl_adv_fct') 1993 ! 1994 ! 1995 ! 1. Initializations 1996 ! ------------------ 1997 ! 1998 IF( neuler == 0 .AND. kt == nit000 ) THEN 1999 z2dtt = rdt 2000 ELSE 2001 z2dtt = 2.0_wp * rdt 2002 ENDIF 2003 ! 2004 zwi(:,:,:) = 0.e0 2005 zwx(:,:,:) = 0.e0 2006 zwy(:,:,:) = 0.e0 2007 ! 2008 ! 2009 ! 2. upstream advection with initial mass fluxes & intermediate update 2010 ! -------------------------------------------------------------------- 2011 IF ( ll_shorizd ) THEN 2012 DO jk = 1, jpkm1 2013 DO jj = 1, jpjm1 2014 DO ji = 1, fs_jpim1 ! vector opt. 2015 ! 2016 zfp_hi = MAX(hu_b(ji,jj) - gdepw_b(ji ,jj ,jk), 0._wp) 2017 zfp_hi = MIN(zfp_hi, e3t_b(ji ,jj ,jk)) 2018 zfp_hi = 0.5_wp *(zfp_hi + SIGN(zfp_hi, zfp_hi-hsmall) ) 2019 ! 2020 zfm_hi = MAX(hu_b(ji,jj) - gdepw_b(ji+1,jj ,jk), 0._wp) 2021 zfm_hi = MIN(zfm_hi, e3t_b(ji+1,jj ,jk)) 2022 zfm_hi = 0.5_wp *(zfm_hi + SIGN(zfm_hi, zfm_hi-hsmall) ) 2023 ! 2024 zfp_hj = MAX(hv_b(ji,jj) - gdepw_b(ji ,jj ,jk), 0._wp) 2025 zfp_hj = MIN(zfp_hj, e3t_b(ji ,jj ,jk)) 2026 zfp_hj = 0.5_wp *(zfp_hj + SIGN(zfp_hj, zfp_hj-hsmall) ) 2027 ! 2028 zfm_hj = MAX(hv_b(ji,jj) - gdepw_b(ji ,jj+1,jk), 0._wp) 2029 zfm_hj = MIN(zfm_hj, e3t_b(ji ,jj+1,jk)) 2030 zfm_hj = 0.5_wp *(zfm_hj + SIGN(zfm_hj, zfm_hj-hsmall) ) 2031 ! 2032 zfp_ui = uin(ji,jj,jk) + ABS( uin(ji,jj,jk) ) 2033 zfm_ui = uin(ji,jj,jk) - ABS( uin(ji,jj,jk) ) 2034 zfp_vj = vin(ji,jj,jk) + ABS( vin(ji,jj,jk) ) 2035 zfm_vj = vin(ji,jj,jk) - ABS( vin(ji,jj,jk) ) 2036 zwx(ji,jj,jk) = 0.5 * e2u(ji,jj) * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk) 2037 zwy(ji,jj,jk) = 0.5 * e1v(ji,jj) * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk) 2038 END DO 2039 END DO 2040 END DO 2041 ELSE 2042 DO jk = 1, jpkm1 2043 DO jj = 1, jpjm1 2044 DO ji = 1, fs_jpim1 ! vector opt. 2045 ! 2046 zfp_hi = e3t_b(ji ,jj ,jk) 2047 zfm_hi = e3t_b(ji+1,jj ,jk) 2048 zfp_hj = e3t_b(ji ,jj ,jk) 2049 zfm_hj = e3t_b(ji ,jj+1,jk) 2050 ! 2051 zfp_ui = uin(ji,jj,jk) + ABS( uin(ji,jj,jk) ) 2052 zfm_ui = uin(ji,jj,jk) - ABS( uin(ji,jj,jk) ) 2053 zfp_vj = vin(ji,jj,jk) + ABS( vin(ji,jj,jk) ) 2054 zfm_vj = vin(ji,jj,jk) - ABS( vin(ji,jj,jk) ) 2055 zwx(ji,jj,jk) = 0.5 * e2u(ji,jj) * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk) 2056 zwy(ji,jj,jk) = 0.5 * e1v(ji,jj) * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk) 2057 END DO 2058 END DO 2059 END DO 2060 ENDIF 2061 2062 ! total advective trend 2063 DO jk = 1, jpkm1 2064 DO jj = 2, jpjm1 2065 DO ji = fs_2, fs_jpim1 ! vector opt. 2066 zbtr = r1_e1e2t(ji,jj) 2067 ! total intermediate advective trends 2068 ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 2069 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) 2070 ! 2071 ! update and guess with monotonic sheme 2072 pta(ji,jj,jk) = pta(ji,jj,jk) + ztra 2073 zwi(ji,jj,jk) = (e3t_b(ji,jj,jk) + z2dtt * ztra ) * tmask(ji,jj,jk) 2074 END DO 2075 END DO 2076 END DO 2077 2078 CALL lbc_lnk( zwi, 'T', 1. ) 2079 2080 IF ( ln_bdy ) THEN 2081 DO ib_bdy=1, nb_bdy 2082 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(1) 2083 ji = idx_bdy(ib_bdy)%nbi(ib,1) 2084 jj = idx_bdy(ib_bdy)%nbj(ib,1) 2085 DO jk = 1, jpkm1 2086 zwi(ji,jj,jk) = e3t_a(ji,jj,jk) 2087 END DO 2088 END DO 2089 END DO 2090 ENDIF 2091 2092 IF ( ln_vvl_dbg ) THEN 2093 zmin = MINVAL( zwi(:,:,:), mask = tmask(:,:,:) == 1.e0 ) 2094 IF( lk_mpp ) CALL mpp_min( zmin ) 2095 IF( zmin < 0._wp) THEN 2096 IF(lwp) CALL ctl_warn('vvl_adv: CFL issue here') 2097 IF(lwp) WRITE(numout,*) zmin 2098 ENDIF 2099 ENDIF 2100 2101 ! 3. antidiffusive flux : high order minus low order 2102 ! -------------------------------------------------- 2103 ! antidiffusive flux on i and j 2104 DO jk = 1, jpkm1 2105 DO jj = 1, jpjm1 2106 DO ji = 1, fs_jpim1 ! vector opt. 2107 zwx(ji,jj,jk) = (e2u(ji,jj) * uin(ji,jj,jk) * e3u_n(ji,jj,jk) & 2108 & - zwx(ji,jj,jk)) * umask(ji,jj,jk) 2109 zwy(ji,jj,jk) = (e1v(ji,jj) * vin(ji,jj,jk) * e3v_n(ji,jj,jk) & 2110 & - zwy(ji,jj,jk)) * vmask(ji,jj,jk) 2111 ! 2112 ! Update advective fluxes 2113 un_td(ji,jj,jk) = un_td(ji,jj,jk) - zwx(ji,jj,jk) 2114 vn_td(ji,jj,jk) = vn_td(ji,jj,jk) - zwy(ji,jj,jk) 2115 END DO 2116 END DO 2117 END DO 2118 2119 CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) ! Lateral boundary conditions 2120 2121 ! 4. monotonicity algorithm 2122 ! ------------------------- 2123 CALL nonosc_2d( e3t_b(:,:,:), zwx, zwy, zwi, z2dtt ) 2124 2125 ! 5. final trend with corrected fluxes 2126 ! ------------------------------------ 2127 ! 2128 ! Update advective fluxes 2129 un_td(:,:,:) = (un_td(:,:,:) + zwx(:,:,:))*umask(:,:,:) 2130 vn_td(:,:,:) = (vn_td(:,:,:) + zwy(:,:,:))*vmask(:,:,:) 2131 ! 2132 DO jk = 1, jpkm1 2133 DO jj = 2, jpjm1 2134 DO ji = fs_2, fs_jpim1 ! vector opt. 2135 ! 2136 zbtr = r1_e1e2t(ji,jj) 2137 ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 2138 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) 2139 ! add them to the general tracer trends 2140 pta(ji,jj,jk) = pta(ji,jj,jk) + ztra 2141 END DO 2142 END DO 2143 END DO 2144 ! 2145 IF( ln_timing ) CALL timing_stop('dom_vvl_adv_fct') 2146 ! 2147 END SUBROUTINE dom_vvl_adv_fct 2148 2149 SUBROUTINE dom_vvl_ups_cor( kt, pta, uin, vin ) 2150 !!---------------------------------------------------------------------- 2151 !! *** ROUTINE dom_vvl_adv_fct *** 2152 !! 2153 !! ** Purpose : Correct for addionnal barotropic fluxes 2154 !! in the upstream direction 2155 !! 2156 !! ** Method : 2157 !! 2158 !! ** Action : - Update diffusive fluxes uin, vin 2159 !! - Remove divergence from thickness tendency 2160 !!---------------------------------------------------------------------- 2161 INTEGER, INTENT( in ) :: kt ! ocean time-step index 2162 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pta ! thickness baroclinic trend 2163 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: uin, vin ! input fluxes 2164 INTEGER :: ji, jj, jk ! dummy loop indices 2165 INTEGER :: ikbu, ikbv, ibot 2166 REAL(wp) :: zbtr, ztra ! local scalar 2167 REAL(wp) :: zdi, zdj ! - - 2168 REAL(wp) :: zfp_hi, zfp_hj ! - - 2169 REAL(wp) :: zfm_hi, zfm_hj ! - - 2170 REAL(wp) :: zfp_ui, zfp_vj ! - - 2171 REAL(wp) :: zfm_ui, zfm_vj ! - - 2172 REAL(wp), DIMENSION(jpi,jpj) :: zbu, zbv, zhu_b, zhv_b 2173 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwx, zwy 2174 !!---------------------------------------------------------------------- 2175 ! 2176 IF( ln_timing ) CALL timing_start('dom_vvl_ups_cor') 2177 ! 2178 ! Compute barotropic flux difference: 2179 zbu(:,:) = 0.e0 2180 zbv(:,:) = 0.e0 2181 DO jj = 1, jpj 2182 DO ji = 1, jpi ! vector opt. 2183 DO jk = 1, jpkm1 2184 zbu(ji,jj) = zbu(ji,jj) - uin(ji,jj,jk) * umask(ji,jj,jk) 2185 zbv(ji,jj) = zbv(ji,jj) - vin(ji,jj,jk) * vmask(ji,jj,jk) 2186 END DO 2187 END DO 2188 ENDDO 2189 2190 ! Compute upstream depths: 2191 zhu_b(:,:) = 0.e0 2192 zhv_b(:,:) = 0.e0 2193 2194 IF ( ll_shorizd ) THEN 2195 ! Correct bottom value 2196 ! considering "shelf horizon depth" 2197 DO jj = 1, jpjm1 2198 DO ji = 1, jpim1 ! vector opt. 2199 zdi = 0.5_wp + 0.5_wp * SIGN(1._wp, zbu(ji,jj)) 2200 zdj = 0.5_wp + 0.5_wp * SIGN(1._wp, zbv(ji,jj)) 2201 DO jk=1, jpkm1 2202 zfp_hi = MAX(hu_b(ji,jj) - gdepw_b(ji ,jj ,jk), 0._wp) 2203 zfp_hi = MIN(zfp_hi, e3t_b(ji ,jj ,jk)) 2204 zfp_hi = 0.5_wp *(zfp_hi + SIGN(zfp_hi, zfp_hi-hsmall) ) 2205 ! 2206 zfm_hi = MAX(hu_b(ji,jj) - gdepw_b(ji+1,jj ,jk), 0._wp) 2207 zfm_hi = MIN(zfm_hi, e3t_b(ji+1,jj ,jk)) 2208 zfm_hi = 0.5_wp *(zfm_hi + SIGN(zfm_hi, zfm_hi-hsmall) ) 2209 ! 2210 zfp_hj = MAX(hv_b(ji,jj) - gdepw_b(ji ,jj ,jk), 0._wp) 2211 zfp_hj = MIN(zfp_hj, e3t_b(ji ,jj ,jk)) 2212 zfp_hj = 0.5_wp *(zfp_hj + SIGN(zfp_hj, zfp_hj-hsmall) ) 2213 ! 2214 zfm_hj = MAX(hv_b(ji,jj) - gdepw_b(ji ,jj+1,jk), 0._wp) 2215 zfm_hj = MIN(zfm_hj, e3t_b(ji ,jj+1,jk)) 2216 zfm_hj = 0.5_wp *(zfm_hj + SIGN(zfm_hj, zfm_hj-hsmall) ) 2217 ! 2218 zhu_b(ji,jj) = zhu_b(ji,jj) + ( zdi * zfp_hi & 2219 & + (1._wp-zdi) * zfm_hi & 2220 & ) * umask(ji,jj,jk) 2221 zhv_b(ji,jj) = zhv_b(ji,jj) + ( zdj * zfp_hj & 2222 & + (1._wp-zdj) * zfm_hj & 2223 & ) * vmask(ji,jj,jk) 2224 END DO 2225 END DO 2226 END DO 2227 ELSE 2228 DO jj = 1, jpjm1 2229 DO ji = 1, jpim1 ! vector opt. 2230 zdi = 0.5_wp + 0.5_wp * SIGN(1._wp, zbu(ji,jj)) 2231 zdj = 0.5_wp + 0.5_wp * SIGN(1._wp, zbv(ji,jj)) 2232 DO jk = 1, jpkm1 2233 zfp_hi = e3t_b(ji ,jj ,jk) 2234 zfm_hi = e3t_b(ji+1,jj ,jk) 2235 zfp_hj = e3t_b(ji ,jj ,jk) 2236 zfm_hj = e3t_b(ji ,jj+1,jk) 2237 ! 2238 zhu_b(ji,jj) = zhu_b(ji,jj) + ( zdi * zfp_hi & 2239 & + (1._wp-zdi) * zfm_hi & 2240 & ) * umask(ji,jj,jk) 2241 ! 2242 zhv_b(ji,jj) = zhv_b(ji,jj) + ( zdj * zfp_hj & 2243 & + (1._wp-zdj) * zfm_hj & 2244 & ) * vmask(ji,jj,jk) 2245 END DO 2246 END DO 2247 END DO 2248 ENDIF 2249 2250 ! Corrective barotropic velocity (times hor. scale factor) 2251 zbu(:,:) = zbu(:,:)/ (zhu_b(:,:)*umask(:,:,1)+1._wp-umask(:,:,1)) 2252 zbv(:,:) = zbv(:,:)/ (zhv_b(:,:)*vmask(:,:,1)+1._wp-vmask(:,:,1)) 2253 2254 CALL lbc_lnk( zbu(:,:), 'U', -1. ) 2255 CALL lbc_lnk( zbv(:,:), 'V', -1. ) 2256 2257 ! Set corrective fluxes in upstream direction: 2258 ! 2259 zwx(:,:,:) = 0.e0 2260 zwy(:,:,:) = 0.e0 2261 IF ( ll_shorizd ) THEN 2262 DO jj = 1, jpjm1 2263 DO ji = 1, fs_jpim1 ! vector opt. 2264 ! upstream scheme 2265 zfp_ui = zbu(ji,jj) + ABS( zbu(ji,jj) ) 2266 zfm_ui = zbu(ji,jj) - ABS( zbu(ji,jj) ) 2267 zfp_vj = zbv(ji,jj) + ABS( zbv(ji,jj) ) 2268 zfm_vj = zbv(ji,jj) - ABS( zbv(ji,jj) ) 2269 DO jk = 1, jpkm1 2270 zfp_hi = MAX(hu_b(ji,jj) - gdepw_b(ji ,jj ,jk), 0._wp) 2271 zfp_hi = MIN(e3t_b(ji ,jj ,jk), zfp_hi) 2272 zfp_hi = 0.5_wp *(zfp_hi + SIGN(zfp_hi, zfp_hi-hsmall) ) 2273 ! 2274 zfm_hi = MAX(hu_b(ji,jj) - gdepw_b(ji+1,jj ,jk), 0._wp) 2275 zfm_hi = MIN(e3t_b(ji+1,jj ,jk), zfm_hi) 2276 zfm_hi = 0.5_wp *(zfm_hi + SIGN(zfm_hi, zfm_hi-hsmall) ) 2277 ! 2278 zfp_hj = MAX(hv_b(ji,jj) - gdepw_b(ji ,jj ,jk), 0._wp) 2279 zfp_hj = MIN(e3t_b(ji ,jj ,jk), zfp_hj) 2280 zfp_hj = 0.5_wp *(zfp_hj + SIGN(zfp_hj, zfp_hj-hsmall) ) 2281 ! 2282 zfm_hj = MAX(hv_b(ji,jj) - gdepw_b(ji ,jj+1,jk), 0._wp) 2283 zfm_hj = MIN(e3t_b(ji ,jj+1,jk), zfm_hj) 2284 zfm_hj = 0.5_wp *(zfm_hj + SIGN(zfm_hj, zfm_hj-hsmall) ) 2285 ! 2286 zwx(ji,jj,jk) = 0.5 * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk) 2287 zwy(ji,jj,jk) = 0.5 * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk) 2288 END DO 2289 END DO 2290 END DO 2291 ELSE 2292 DO jj = 1, jpjm1 2293 DO ji = 1, fs_jpim1 ! vector opt. 2294 ! upstream scheme 2295 zfp_ui = zbu(ji,jj) + ABS( zbu(ji,jj) ) 2296 zfm_ui = zbu(ji,jj) - ABS( zbu(ji,jj) ) 2297 zfp_vj = zbv(ji,jj) + ABS( zbv(ji,jj) ) 2298 zfm_vj = zbv(ji,jj) - ABS( zbv(ji,jj) ) 2299 DO jk = 1, jpkm1 2300 zfp_hi = e3t_b(ji ,jj ,jk) 2301 zfm_hi = e3t_b(ji+1,jj ,jk) 2302 zfp_hj = e3t_b(ji ,jj ,jk) 2303 zfm_hj = e3t_b(ji ,jj+1,jk) 2304 ! 2305 zwx(ji,jj,jk) = 0.5 * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk) 2306 zwy(ji,jj,jk) = 0.5 * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk) 2307 END DO 2308 END DO 2309 END DO 2310 ENDIF 2311 CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) ! Lateral boundary conditions 2312 2313 uin(:,:,:) = uin(:,:,:) + zwx(:,:,:) 2314 vin(:,:,:) = vin(:,:,:) + zwy(:,:,:) 2315 ! 2316 ! Update trend with corrective fluxes: 2317 DO jk = 1, jpkm1 2318 DO jj = 2, jpjm1 2319 DO ji = fs_2, fs_jpim1 ! vector opt. 2320 ! 2321 zbtr = r1_e1e2t(ji,jj) 2322 ! total advective trends 2323 ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 2324 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) 2325 ! add them to the general tracer trends 2326 pta(ji,jj,jk) = pta(ji,jj,jk) + ztra 2327 END DO 2328 END DO 2329 END DO 2330 ! 2331 IF( ln_timing ) CALL timing_stop('dom_vvl_ups_cor') 2332 ! 2333 END SUBROUTINE dom_vvl_ups_cor 2334 2335 SUBROUTINE nonosc_2d( pbef, paa, pbb, paft, p2dt ) 2336 !!--------------------------------------------------------------------- 2337 !! *** ROUTINE nonosc_2d *** 2338 !! 2339 !! ** Purpose : compute monotonic thickness fluxes from the upstream 2340 !! scheme and the before field by a nonoscillatory algorithm 2341 !! 2342 !! ** Method : ... ??? 2343 !! warning : pbef and paft must be masked, but the boundaries 2344 !! conditions on the fluxes are not necessary zalezak (1979) 2345 !! drange (1995) multi-dimensional forward-in-time and upstream- 2346 !! in-space based differencing for fluid 2347 !!---------------------------------------------------------------------- 2348 ! 2349 !!---------------------------------------------------------------------- 2350 REAL(wp) , INTENT(in ) :: p2dt ! vertical profile of tracer time-step 2351 REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in ) :: pbef, paft ! before & after field 2352 REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) :: paa, pbb ! monotonic fluxes in the 3 directions 2353 ! 2354 INTEGER :: ji, jj, jk ! dummy loop indices 2355 REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt ! local scalars 2356 REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo ! - - 2357 REAL(wp) :: zupip1, zupim1, zupjp1, zupjm1, zupb, zupa 2358 REAL(wp) :: zdoip1, zdoim1, zdojp1, zdojm1, zdob, zdoa 2359 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zbetup, zbetdo, zbup, zbdo 2360 !!---------------------------------------------------------------------- 2361 ! 2362 IF( ln_timing ) CALL timing_start('nonosc2') 2363 ! 2364 zbig = 1.e+40_wp 2365 zrtrn = 1.e-15_wp 2366 zbetup(:,:,jpk) = 0._wp ; zbetdo(:,:,jpk) = 0._wp 2367 2368 2369 ! Search local extrema 2370 ! -------------------- 2371 ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land 2372 zbup = MAX( pbef * tmask - zbig * ( 1.e0 - tmask ), & 2373 & paft * tmask - zbig * ( 1.e0 - tmask ) ) 2374 zbdo = MIN( pbef * tmask + zbig * ( 1.e0 - tmask ), & 2375 & paft * tmask + zbig * ( 1.e0 - tmask ) ) 2376 2377 DO jk = 1, jpkm1 2378 z2dtt = p2dt 2379 DO jj = 2, jpjm1 2380 DO ji = fs_2, fs_jpim1 ! vector opt. 2381 2382 ! search maximum in neighbourhood 2383 zup = MAX( zbup(ji ,jj ,jk ), & 2384 & zbup(ji-1,jj ,jk ), zbup(ji+1,jj ,jk ), & 2385 & zbup(ji ,jj-1,jk ), zbup(ji ,jj+1,jk )) 2386 2387 ! search minimum in neighbourhood 2388 zdo = MIN( zbdo(ji ,jj ,jk ), & 2389 & zbdo(ji-1,jj ,jk ), zbdo(ji+1,jj ,jk ), & 2390 & zbdo(ji ,jj-1,jk ), zbdo(ji ,jj+1,jk )) 2391 2392 ! positive part of the flux 2393 zpos = MAX( 0., paa(ji-1,jj ,jk ) ) - MIN( 0., paa(ji ,jj ,jk ) ) & 2394 & + MAX( 0., pbb(ji ,jj-1,jk ) ) - MIN( 0., pbb(ji ,jj ,jk ) ) 2395 2396 ! negative part of the flux 2397 zneg = MAX( 0., paa(ji ,jj ,jk ) ) - MIN( 0., paa(ji-1,jj ,jk ) ) & 2398 & + MAX( 0., pbb(ji ,jj ,jk ) ) - MIN( 0., pbb(ji ,jj-1,jk ) ) 2399 2400 ! up & down beta terms 2401 zbt = e1t(ji,jj) * e2t(ji,jj) / z2dtt 2402 zbetup(ji,jj,jk) = ( zup - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt 2403 zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo ) / ( zneg + zrtrn ) * zbt 2404 END DO 2405 END DO 2406 END DO 2407 2408 CALL lbc_lnk( zbetup, 'T', 1. ) ; CALL lbc_lnk( zbetdo, 'T', 1. ) ! lateral boundary cond. (unchanged sign) 2409 2410 ! 3. monotonic flux in the i & j direction (paa & pbb) 2411 ! ---------------------------------------- 2412 DO jk = 1, jpkm1 2413 DO jj = 2, jpjm1 2414 DO ji = fs_2, fs_jpim1 ! vector opt. 2415 zau = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) ) 2416 zbu = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) ) 2417 zcu = ( 0.5 + SIGN( 0.5 , paa(ji,jj,jk) ) ) 2418 paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1.e0 - zcu) * zbu ) 2419 2420 zav = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) ) 2421 zbv = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) ) 2422 zcv = ( 0.5 + SIGN( 0.5 , pbb(ji,jj,jk) ) ) 2423 pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1.e0 - zcv) * zbv ) 2424 END DO 2425 END DO 2426 END DO 2427 CALL lbc_lnk( paa, 'U', -1. ) ; CALL lbc_lnk( pbb, 'V', -1. ) ! lateral boundary condition (changed sign) 2428 ! 2429 IF( ln_timing ) CALL timing_stop('nonosc2') 2430 ! 2431 END SUBROUTINE nonosc_2d 2432 1245 2433 !!====================================================================== 1246 2434 END MODULE domvvl -
NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/DYN/dynnxt.F90
r9598 r10116 95 95 ! 96 96 INTEGER :: ji, jj, jk ! dummy loop indices 97 INTEGER :: ikt 97 INTEGER :: ikt, jkbot ! local integers 98 98 REAL(wp) :: zue3a, zue3n, zue3b, zuf, zcoef ! local scalars 99 99 REAL(wp) :: zve3a, zve3n, zve3b, zvf, z1_2dt ! - - 100 REAL(wp) :: zhcri, zmsku, zwgtu, zmskv, zwgtv ! - - 100 101 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zue, zve 101 102 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ze3u_f, ze3v_f, zua, zva … … 136 137 END DO 137 138 ENDIF 139 ENDIF 140 141 IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN 142 ! ! Duplicate baroclinic velocities from above in massless layers 143 ! ! Correction to handle correct barotropic velocities right after 144 zhcri = 0.1_wp 145 DO ji = 2, jpim1 146 DO jj= 2, jpjm1 147 jkbot=jpkm1 148 DO jk=jpkm1,2,-1 149 IF (e3u_a(ji,jj,jk)>zhcri) jkbot = jk+1 150 END DO 151 DO jk=jkbot,jpkm1 152 zwgtu = MIN(zhcri, e3u_a(ji,jj,jk)) 153 zmsku = 0.5_wp * ( 1._wp + SIGN(1._wp, e3u_a(ji,jj,jk) - 0.01_wp) ) 154 ua(ji,jj,jk) = zmsku * (zwgtu * ua(ji,jj,jk) + (zhcri - zwgtu) * ua(ji,jj,jk-1)*zwgtu/(zwgtu+0.01_wp)) / zhcri 155 END DO 156 jkbot=jpkm1 157 DO jk=jpkm1,2,-1 158 IF (e3v_a(ji,jj,jk)>zhcri) jkbot = jk+1 159 END DO 160 DO jk=jkbot,jpkm1 161 zwgtv = MIN(zhcri, e3v_a(ji,jj,jk)) 162 zmskv = 0.5_wp * ( 1._wp + SIGN(1._wp, e3v_a(ji,jj,jk) - 0.01_wp) ) 163 va(ji,jj,jk) = zmskv * (zwgtv * va(ji,jj,jk) + (zhcri - zwgtv) * va(ji,jj,jk-1)*zwgtv/(zwgtv+0.01_wp)) / zhcri 164 END DO 165 END DO 166 END DO 138 167 ENDIF 139 168 … … 291 320 zuf = ( zue3n + atfp * ( zue3b - 2._wp * zue3n + zue3a ) ) / ze3u_f(ji,jj,jk) 292 321 zvf = ( zve3n + atfp * ( zve3b - 2._wp * zve3n + zve3a ) ) / ze3v_f(ji,jj,jk) 322 zmsku = 0.5_wp * ( 1._wp + SIGN(1._wp, ze3u_f(ji,jj,jk) - 0.01_wp) ) 323 zmskv = 0.5_wp * ( 1._wp + SIGN(1._wp, ze3v_f(ji,jj,jk) - 0.01_wp) ) 293 324 ! 294 ub(ji,jj,jk) = z uf! ub <-- filtered velocity295 vb(ji,jj,jk) = z vf325 ub(ji,jj,jk) = zmsku * zuf ! ub <-- filtered velocity 326 vb(ji,jj,jk) = zmskv * zvf 296 327 un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua 297 328 vn(ji,jj,jk) = va(ji,jj,jk) -
NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/src/OCE/SBC/sbcwave.F90
r9966 r10116 21 21 USE zdf_oce, ONLY : ln_zdfswm 22 22 USE bdy_oce ! open boundary condition variables 23 USE domvvl ! domain: variable volume layers24 23 ! 25 24 USE iom ! I/O manager library
Note: See TracChangeset
for help on using the changeset viewer.
