- Timestamp:
- 2018-02-23T16:52:00+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2018/dev_r8864_nemo_v3_6_ZTILDE/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r6348 r9353 9 9 !! vvl option includes z_star and z_tilde coordinates 10 10 !! 3.6 ! 2014-11 (P. Mathiot) add ice shelf capability 11 !! ! 2018-01 (J. Chanut) improve ztilde robustness 11 12 !!---------------------------------------------------------------------- 12 13 !! 'key_vvl' variable volume … … 33 34 USE wrk_nemo ! Memory allocation 34 35 USE timing ! Timing 36 USE bdy_oce ! ocean open boundary conditions 35 37 36 38 IMPLICIT NONE … … 48 50 LOGICAL , PUBLIC :: ln_vvl_ztilde = .FALSE. ! ztilde vertical coordinate 49 51 LOGICAL , PUBLIC :: ln_vvl_layer = .FALSE. ! level vertical coordinate 50 LOGICAL , PUBLIC :: ln_vvl_ztilde_as_zstar = .FALSE. ! ztilde vertical coordinate 51 LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor = .FALSE. ! ztilde vertical coordinate 52 LOGICAL , PUBLIC :: ln_vvl_kepe = .FALSE. ! kinetic/potential energy transfer 52 LOGICAL :: ln_vvl_ztilde_as_zstar = .FALSE. ! ztilde vertical coordinate 53 LOGICAL :: ln_vvl_zstar_at_eqtor = .FALSE. ! revert to zstar at equator 54 LOGICAL :: ln_vvl_zstar_on_shelf =.FALSE. ! revert to zstar on shelves 55 LOGICAL :: ln_vvl_kepe =.FALSE. ! kinetic/potential energy transfer 56 LOGICAL :: ln_vvl_adv_fct =.FALSE. ! Centred thickness advection 57 LOGICAL :: ln_vvl_adv_cn2 =.TRUE. ! FCT thickness advection 58 LOGICAL :: ln_vvl_dbg = .FALSE. ! debug control prints 59 LOGICAL :: ln_vvl_lap ! Laplacian thickness diffusion 60 LOGICAL :: ln_vvl_blp ! Bilaplacian thickness diffusion 61 LOGICAL :: ln_vvl_regrid ! ensure layer separation 53 62 ! ! conservation: not used yet 54 REAL(wp) :: rn_ahe3 ! thickness diffusion coefficient 63 REAL(wp) :: rn_ahe3_lap ! thickness diffusion coefficient (Laplacian) 64 REAL(wp) :: rn_ahe3_blp ! thickness diffusion coefficient (Bilaplacian) 55 65 REAL(wp) :: rn_rst_e3t ! ztilde to zstar restoration timescale [days] 56 66 REAL(wp) :: rn_lf_cutoff ! cutoff frequency for low-pass filter [days] 57 67 REAL(wp) :: rn_zdef_max ! maximum fractional e3t deformation 58 LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE. ! debug control prints59 68 60 69 !! * Module variables 61 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td ! thickness diffusion transport 62 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf ! low frequency part of hz divergence 63 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n ! baroclinic scale factors 64 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a ! baroclinic scale factors 65 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_e3t ! retoring period for scale factors 66 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_hdv ! retoring period for low freq. divergence 70 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td ! thickness diffusion transport 71 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdivn_lf, un_lf, vn_lf ! 1st order filtered arrays 72 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n ! baroclinic scale factors 73 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a, dtilde_e3t_a ! baroclinic scale factors 74 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_e3t ! restoring period for scale factors 75 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_hdv ! restoring period for low freq. divergence 76 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: tildemask ! mask tilde tendency 77 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: hsm, dsm ! 78 INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:) :: i_int_bot 67 79 68 80 !! * Substitutions … … 85 97 ALLOCATE( tilde_e3t_b(jpi,jpj,jpk) , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) , & 86 98 & dtilde_e3t_a(jpi,jpj,jpk) , un_td (jpi,jpj,jpk) , vn_td (jpi,jpj,jpk) , & 99 & i_int_bot(jpi,jpj), tildemask(jpi,jpj), hsm(jpi,jpj), dsm(jpi,jpj), & 87 100 & STAT = dom_vvl_alloc ) 88 101 IF( lk_mpp ) CALL mpp_sum ( dom_vvl_alloc ) … … 92 105 ENDIF 93 106 IF( ln_vvl_ztilde ) THEN 94 ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdiv_lf(jpi,jpj,jpk) , STAT= dom_vvl_alloc ) 107 ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj), & 108 & hdivn_lf(jpi,jpj,jpk), & 109 & un_lf(jpi,jpj,jpk), & 110 & vn_lf(jpi,jpj,jpk), STAT= dom_vvl_alloc ) 95 111 IF( lk_mpp ) CALL mpp_sum ( dom_vvl_alloc ) 96 112 IF( dom_vvl_alloc /= 0 ) CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') … … 124 140 USE phycst, ONLY : rpi, rsmall, rad 125 141 !! * Local declarations 126 INTEGER :: ji, jj,jk142 INTEGER :: ji, jj, jk 127 143 INTEGER :: ii0, ii1, ij0, ij1 128 REAL(wp):: zcoef 144 REAL(wp):: zcoef, zwgt, ztmp, zhmin, zhmax 129 145 !!---------------------------------------------------------------------- 130 146 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_init') … … 155 171 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' ) 156 172 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' ) 157 CALL dom_vvl_interpol( fse3 u_n(:,:,:), fse3f_n(:,:,:), 'F' )173 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3f_n(:,:,:), 'F' ) 158 174 ! Vertical scale factor interpolations 159 175 ! ------------------------------------ … … 199 215 hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 200 216 END DO 201 hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1. - umask_i(:,:) )202 hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1. - vmask_i(:,:) )217 hur_b(:,:) = umask_i(:,:) / ( hu_b(:,:) + 1._wp - umask_i(:,:) ) 218 hvr_b(:,:) = vmask_i(:,:) / ( hv_b(:,:) + 1._wp - vmask_i(:,:) ) 203 219 204 220 ! Restoring frequencies for z_tilde coordinate 205 221 ! ============================================ 222 tildemask(:,:) = 1._wp 223 206 224 IF( ln_vvl_ztilde ) THEN 207 225 ! Values in days provided via the namelist; use rsmall to avoid possible division by zero errors with faulty settings 208 226 frq_rst_e3t(:,:) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.0_wp ) 209 227 frq_rst_hdv(:,:) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp ) 228 ! tendency mask: 229 ! 210 230 IF( ln_vvl_ztilde_as_zstar ) THEN 211 231 ! Ignore namelist settings and use these next two to emulate z-star using z-tilde 212 232 frq_rst_e3t(:,:) = 0.0_wp 213 233 frq_rst_hdv(:,:) = 1.0_wp / rdt 234 tildemask(:,:) = 0._wp 214 235 ENDIF 236 215 237 IF ( ln_vvl_zstar_at_eqtor ) THEN 216 238 DO jj = 1, jpj … … 219 241 ! values outside the equatorial band and transition zone (ztilde) 220 242 frq_rst_e3t(ji,jj) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.e0_wp ) 221 frq_rst_hdv(ji,jj) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 243 ! frq_rst_hdv(ji,jj) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 244 222 245 ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN 223 246 ! values inside the equatorial band (ztilde as zstar) 224 247 frq_rst_e3t(ji,jj) = 0.0_wp 225 frq_rst_hdv(ji,jj) = 1.0_wp / rdt 248 ! frq_rst_hdv(ji,jj) = 1.0_wp / rdt 249 tildemask(ji,jj) = 0._wp 226 250 ELSE 227 251 ! values in the transition band (linearly vary from ztilde to ztilde as zstar values) … … 229 253 & * ( 1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & 230 254 & * 180._wp / 3.5_wp ) ) 231 frq_rst_hdv(ji,jj) = (1.0_wp / rdt) & 232 & + ( frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp & 233 & * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & 234 & * 180._wp / 3.5_wp ) ) 255 ! frq_rst_hdv(ji,jj) = (1.0_wp / rdt) & 256 ! & + ( frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp & 257 ! & * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & 258 ! & * 180._wp / 3.5_wp ) ) 259 tildemask(ji,jj) = 0.5_wp * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & 260 & * 180._wp / 3.5_wp ) ) 235 261 ENDIF 236 262 END DO … … 240 266 ij0 = 128 ; ij1 = 135 ; 241 267 frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.0_wp 242 frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0_wp / rdt 268 tildemask( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.0_wp 269 ! frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0_wp / rdt 243 270 ENDIF 244 271 ENDIF 272 ! 273 IF ( ln_vvl_zstar_on_shelf ) THEN 274 zhmin = 50._wp 275 zhmax = 100._wp 276 DO jj = 1, jpj 277 DO ji = 1, jpi 278 zwgt = 1._wp 279 IF(( ht_0(ji,jj)>zhmin).AND.(ht_0(ji,jj) <=zhmax)) THEN 280 zwgt = (ht_0(ji,jj)-zhmin)/(zhmax-zhmin) 281 ELSEIF ( ht_0(ji,jj)<=zhmin) THEN 282 zwgt = 0._wp 283 ENDIF 284 frq_rst_e3t(ji,jj) = MIN(frq_rst_e3t(ji,jj), frq_rst_e3t(ji,jj)*zwgt) 285 tildemask(ji,jj) = MIN(tildemask(ji,jj), zwgt) 286 END DO 287 END DO 288 ENDIF 289 ! 290 ztmp = MAXVAL( frq_rst_hdv(:,:) ) 291 IF( lk_mpp ) CALL mpp_max( ztmp ) ! max over the global domain 292 ! 293 IF ( (ztmp*rdt) > 1._wp) CALL ctl_stop( 'dom_vvl_init: rn_lf_cuttoff is too small' ) 294 ! 245 295 ENDIF 246 296 … … 272 322 !! Reference : Leclair, M., and Madec, G. 2011, Ocean Modelling. 273 323 !!---------------------------------------------------------------------- 274 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t 324 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t, ztu, ztv 275 325 REAL(wp), POINTER, DIMENSION(:,: ) :: zht, z_scale, zwu, zwv, zhdiv 276 326 !! * Arguments … … 279 329 !! * Local declarations 280 330 INTEGER :: ji, jj, jk ! dummy loop indices 331 INTEGER :: ib, ib_bdy, ip, jp ! " " " 281 332 INTEGER , DIMENSION(3) :: ijk_max, ijk_min ! temporary integers 282 333 REAL(wp) :: z2dt ! temporary scalars 283 334 REAL(wp) :: z_tmin, z_tmax ! temporary scalars 335 REAL(wp) :: zalpha, zwgt ! temporary scalars 336 REAL(wp) :: zdu, zdv 284 337 LOGICAL :: ll_do_bclinic ! temporary logical 285 338 !!---------------------------------------------------------------------- 286 339 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_nxt') 287 340 CALL wrk_alloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv ) 288 CALL wrk_alloc( jpi, jpj, jpk, ze3t 341 CALL wrk_alloc( jpi, jpj, jpk, ze3t, ztu, ztv ) 289 342 290 343 IF(kt == nit000) THEN … … 296 349 ll_do_bclinic = .TRUE. 297 350 IF( PRESENT(kcall) ) THEN 298 IF ( kcall == 2 .AND. ln_vvl_ztilde ) ll_do_bclinic = .FALSE.351 IF ( kcall == 2 .AND. ln_vvl_ztilde.OR.ln_vvl_layer ) ll_do_bclinic = .FALSE. 299 352 ENDIF 300 353 … … 307 360 ! ! --------------------------------------------- ! 308 361 309 z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) )362 z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1._wp - ssmask(:,:) ) 310 363 DO jk = 1, jpkm1 311 364 ! formally this is the same as fse3t_a = e3t_0*(1+ssha/ht_0) 312 365 fse3t_a(:,:,jk) = fse3t_b(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 313 366 END DO 314 315 IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate ! 316 ! ! ------baroclinic part------ ! 317 318 ! I - initialization 319 ! ================== 320 321 ! 1 - barotropic divergence 322 ! ------------------------- 367 368 IF((ln_vvl_ztilde .OR. ln_vvl_layer).AND.ll_do_bclinic ) THEN ! z_tilde or layer coordinate ! 369 370 tilde_e3t_a(:,:,:) = 0.0_wp ! tilde_e3t_a used to store tendency terms 371 un_td(:,:,:) = 0.0_wp ! Transport corrections 372 vn_td(:,:,:) = 0.0_wp 373 323 374 zhdiv(:,:) = 0. 324 zht(:,:) = 0.325 375 DO jk = 1, jpkm1 326 376 zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk) 327 zht (:,:) = zht (:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 328 END DO 329 zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) ) 330 331 ! 2 - Low frequency baroclinic horizontal divergence (z-tilde case only) 332 ! -------------------------------------------------- 333 IF( ln_vvl_ztilde ) THEN 334 IF( kt .GT. nit000 ) THEN 377 END DO 378 zhdiv(:,:) = zhdiv(:,:) / ( ht_0(:,:) + sshn(:,:) + 1._wp - tmask_i(:,:) ) 379 380 ! Thickness advection: 381 ! -------------------- 382 ! Set advection velocities and source term 383 IF ( ln_vvl_ztilde ) THEN 384 ! 385 IF ((kt==nit000).AND.(neuler==0)) THEN 335 386 DO jk = 1, jpkm1 336 hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:) & 337 & * ( hdiv_lf(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) ) 387 ztu(:,:,jk) = un(:,:,jk) 388 ztv(:,:,jk) = vn(:,:,jk) 389 END DO 390 ELSE 391 DO jk = 1, jpkm1 392 ztu(:,:,jk) = (un(:,:,jk)-un_lf(:,:,jk)/fse3u_n(:,:,jk)*r1_e2u(:,:))*umask(:,:,jk) 393 ztv(:,:,jk) = (vn(:,:,jk)-vn_lf(:,:,jk)/fse3v_n(:,:,jk)*r1_e1v(:,:))*vmask(:,:,jk) 394 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) + hdivn_lf(:,:,jk) 338 395 END DO 339 396 ENDIF 340 END IF 341 342 ! II - after z_tilde increments of vertical scale factors 343 ! ======================================================= 344 tilde_e3t_a(:,:,:) = 0.0_wp ! tilde_e3t_a used to store tendency terms 345 346 ! 1 - High frequency divergence term 347 ! ---------------------------------- 348 IF( ln_vvl_ztilde ) THEN ! z_tilde case 397 ! 398 ELSEIF ( ln_vvl_layer ) THEN 399 ! 349 400 DO jk = 1, jpkm1 350 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 351 END DO 352 ELSE ! layer case 401 ztu(:,:,jk) = un(:,:,jk) 402 ztv(:,:,jk) = vn(:,:,jk) 403 END DO 404 ! 405 ENDIF 406 ! 407 ! Do advection 408 IF (ln_vvl_adv_fct) THEN 409 CALL dom_vvl_adv_fct( kt, tilde_e3t_a, ztu, ztv ) 410 ! 411 ELSEIF (ln_vvl_adv_cn2) THEN 353 412 DO jk = 1, jpkm1 354 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 355 END DO 356 END IF 357 358 ! 2 - Restoring term (z-tilde case only) 359 ! ------------------ 413 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - fse3t_n(:,:,jk) * hdivn(:,:,jk) 414 END DO 415 ENDIF 416 ! 417 ! Thickness anomlaly diffusion: 418 ! ----------------------------- 419 zwu(:,:) = 0.0_wp 420 zwv(:,:) = 0.0_wp 421 ztu(:,:,:) = 0.0_wp 422 ztv(:,:,:) = 0.0_wp 423 424 IF ( ln_vvl_blp ) THEN ! Bilaplacian 425 DO jk = 1, jpkm1 426 DO jj = 1, jpjm1 ! First derivative (gradient) 427 DO ji = 1, fs_jpim1 ! vector opt. 428 ztu(ji,jj,jk) = umask(ji,jj,jk) * re2u_e1u(ji,jj) & 429 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj ,jk) ) 430 ztv(ji,jj,jk) = vmask(ji,jj,jk) * re1v_e2v(ji,jj) & 431 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji ,jj+1,jk) ) 432 END DO 433 END DO 434 435 DO jj = 2, jpjm1 ! Second derivative (divergence) time the eddy diffusivity coefficient 436 DO ji = fs_2, fs_jpim1 ! vector opt. 437 zht(ji,jj) = rn_ahe3_blp * r1_e12t(ji,jj) * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) & 438 & + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) 439 END DO 440 END DO 441 442 #if defined key_bdy 443 DO ib_bdy=1, nb_bdy 444 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(1) 445 ji = idx_bdy(ib_bdy)%nbi(ib,1) 446 jj = idx_bdy(ib_bdy)%nbj(ib,1) 447 zht(ji,jj) = 0._wp 448 END DO 449 END DO 450 #endif 451 CALL lbc_lnk( zht, 'T', 1. ) ! Lateral boundary conditions (unchanged sgn) 452 453 DO jj = 1, jpjm1 ! third derivative (gradient) 454 DO ji = 1, fs_jpim1 ! vector opt. 455 ztu(ji,jj,jk) = umask(ji,jj,jk) * re2u_e1u(ji,jj) * ( zht(ji+1,jj ) - zht(ji,jj) ) 456 ztv(ji,jj,jk) = vmask(ji,jj,jk) * re1v_e2v(ji,jj) * ( zht(ji ,jj+1) - zht(ji,jj) ) 457 zwu(ji,jj) = zwu(ji,jj) + ztu(ji,jj,jk) 458 zwv(ji,jj) = zwv(ji,jj) + ztv(ji,jj,jk) 459 END DO 460 END DO 461 END DO 462 ENDIF 463 464 IF ( ln_vvl_lap ) THEN ! Laplacian 465 DO jk = 1, jpkm1 ! First derivative (gradient) 466 DO jj = 1, jpjm1 467 DO ji = 1, fs_jpim1 ! vector opt. 468 zdu = rn_ahe3_lap * umask(ji,jj,jk) * re2u_e1u(ji,jj) & 469 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj ,jk) ) 470 zdv = rn_ahe3_lap * vmask(ji,jj,jk) * re1v_e2v(ji,jj) & 471 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji ,jj+1,jk) ) 472 zwu(ji,jj) = zwu(ji,jj) + zdu 473 zwv(ji,jj) = zwv(ji,jj) + zdv 474 ztu(ji,jj,jk) = ztu(ji,jj,jk) + zdu 475 ztv(ji,jj,jk) = ztv(ji,jj,jk) + zdv 476 END DO 477 END DO 478 END DO 479 ENDIF 480 481 ! Ensure barotropic fluxes are null: 482 ! DO jj = 1, jpj 483 ! DO ji = 1, jpi 484 ! DO jk = 1, jpkm1 485 ! ztu(ji,jj,jk) = ztu(ji,jj,jk) - zwu(ji,jj)*fse3u_b(ji,jj,jk)*hur_b(ji,jj)*umask(ji,jj,jk) 486 ! ztv(ji,jj,jk) = ztv(ji,jj,jk) - zwv(ji,jj)*fse3v_b(ji,jj,jk)*hvr_b(ji,jj)*vmask(ji,jj,jk) 487 ! END DO 488 ! END DO 489 ! END DO 490 DO jj = 1, jpj 491 DO ji = 1, jpi 492 ztu(ji,jj,mbku(ji,jj)) = ztu(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 493 ztv(ji,jj,mbkv(ji,jj)) = ztv(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 494 END DO 495 END DO 496 497 ! divergence of diffusive fluxes 498 DO jk = 1, jpkm1 499 DO jj = 2, jpjm1 500 DO ji = fs_2, fs_jpim1 ! vector opt. 501 tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + ( ztu(ji-1,jj ,jk) - ztu(ji,jj,jk) & 502 & + ztv(ji ,jj-1,jk) - ztv(ji,jj,jk) & 503 & ) * r1_e12t(ji,jj) 504 END DO 505 END DO 506 END DO 507 508 un_td(:,:,:) = un_td(:,:,:) + ztu(:,:,:) 509 vn_td(:,:,:) = vn_td(:,:,:) + ztv(:,:,:) 510 CALL lbc_lnk( un_td , 'U' , -1.) 511 CALL lbc_lnk( vn_td , 'V' , -1.) 512 ! 513 ! 514 ! Restoring: 360 515 IF( ln_vvl_ztilde ) THEN 361 516 DO jk = 1, jpk 362 517 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk) 363 518 END DO 364 END IF 365 366 ! 3 - Thickness diffusion term 367 ! ---------------------------- 368 zwu(:,:) = 0.0_wp 369 zwv(:,:) = 0.0_wp 370 ! a - first derivative: diffusive fluxes 519 ENDIF 520 521 ! Remove "external thickness" tendency: 371 522 DO jk = 1, jpkm1 372 DO jj = 1, jpjm1 373 DO ji = 1, fs_jpim1 ! vector opt. 374 un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * re2u_e1u(ji,jj) & 375 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj ,jk) ) 376 vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * re1v_e2v(ji,jj) & 377 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji ,jj+1,jk) ) 378 zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 379 zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 380 END DO 381 END DO 382 END DO 383 ! b - correction for last oceanic u-v points 384 DO jj = 1, jpj 385 DO ji = 1, jpi 386 un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 387 vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 388 END DO 389 END DO 390 ! c - second derivative: divergence of diffusive fluxes 391 DO jk = 1, jpkm1 392 DO jj = 2, jpjm1 393 DO ji = fs_2, fs_jpim1 ! vector opt. 394 tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + ( un_td(ji-1,jj ,jk) - un_td(ji,jj,jk) & 395 & + vn_td(ji ,jj-1,jk) - vn_td(ji,jj,jk) & 396 & ) * r1_e12t(ji,jj) 397 END DO 398 END DO 399 END DO 400 ! d - thickness diffusion transport: boundary conditions 401 ! (stored for tracer advction and continuity equation) 402 CALL lbc_lnk( un_td , 'U' , -1._wp) 403 CALL lbc_lnk( vn_td , 'V' , -1._wp) 404 405 ! 4 - Time stepping of baroclinic scale factors 406 ! --------------------------------------------- 523 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) + fse3t_n(:,:,jk) * zhdiv(:,:) 524 END DO 525 407 526 ! Leapfrog time stepping 408 527 ! ~~~~~~~~~~~~~~~~~~~~~~ … … 412 531 z2dt = 2.0_wp * rdt 413 532 ENDIF 414 CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1._wp ) 533 415 534 tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:) 535 536 ! Revert to zstar locally: 537 ! ~~~~~~~~~~~~~~~~~~~~~~~~ 538 DO jk=1,jpkm1 539 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) * tildemask(:,:) 540 END DO 541 542 ! Ensure layer separation: 543 ! ~~~~~~~~~~~~~~~~~~~~~~~~ 544 IF ( ln_vvl_regrid ) CALL dom_vvl_regrid( kt ) 545 546 ! Boundary conditions: 547 ! ~~~~~~~~~~~~~~~~~~~~ 548 #if defined key_bdy 549 DO ib_bdy=1, nb_bdy 550 DO ib = 1, idx_bdy(ib_bdy)%nblenrim(1) 551 !! DO ib = 1, idx_bdy(ib_bdy)%nblen(1) 552 ji = idx_bdy(ib_bdy)%nbi(ib,1) 553 jj = idx_bdy(ib_bdy)%nbj(ib,1) 554 zwgt = idx_bdy(ib_bdy)%nbw(ib,1) 555 ip = bdytmask(ji+1,jj ) - bdytmask(ji-1,jj ) 556 jp = bdytmask(ji ,jj+1) - bdytmask(ji ,jj-1) 557 DO jk = 1, jpkm1 558 tilde_e3t_a(ji,jj,jk) = 0.e0 559 !! tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) * (1._wp - zwgt) 560 !! tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji+ip,jj+jp,jk) * tmask(ji+ip,jj+jp,jk) 561 END DO 562 END DO 563 END DO 564 #endif 565 CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1. ) 416 566 417 567 ! Maximum deformation control … … 446 596 ENDIF 447 597 ENDIF 448 ! - ML - end test 449 ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below 450 tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:), rn_zdef_max * e3t_0(:,:,:) ) 451 tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) ) 452 453 ! 454 ! "tilda" change in the after scale factor 455 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 598 599 IF ( ln_vvl_ztilde ) THEN 600 zalpha = rdt * 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp ) 601 DO jk = 1, jpkm1 602 ztu(:,:,jk) = un(:,:,jk) * fse3u_n(:,:,jk) * e2u(:,:) + un_td(:,:,jk) 603 ztv(:,:,jk) = vn(:,:,jk) * fse3v_n(:,:,jk) * e1v(:,:) + vn_td(:,:,jk) 604 ze3t(:,:,jk) = -fse3t_n(:,:,jk) * zhdiv(:,:) 605 ! 606 ! DO jj = 2, jpjm1 607 ! DO ji = fs_2, fs_jpim1 ! vector opt. 608 ! 609 ! ze3t(ji,jj,jk) = -fse3t_n(ji,jj,jk) * zhdiv(ji,jj) & 610 ! & + ( un_td(ji,jj,jk) - un_td(ji-1,jj ,jk) & 611 ! & + vn_td(ji,jj,jk) - vn_td(ji ,jj-1,jk) ) & 612 ! & / ( e1t(ji,jj) * e2t(ji,jj) ) 613 ! END DO 614 ! END DO 615 END DO 616 ! 617 un_lf(:,:,:) = un_lf(:,:,:) * (1._wp - zalpha) + zalpha * ztu(:,:,:) 618 vn_lf(:,:,:) = vn_lf(:,:,:) * (1._wp - zalpha) + zalpha * ztv(:,:,:) 619 hdivn_lf(:,:,:) = hdivn_lf(:,:,:) * (1._wp - zalpha) + zalpha * ze3t(:,:,:) 620 ENDIF 621 622 ENDIF 623 624 IF((ln_vvl_ztilde .OR. ln_vvl_layer).AND.(.NOT.ll_do_bclinic) ) THEN 625 zhdiv(:,:) = 0. 456 626 DO jk = 1, jpkm1 457 dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk) 458 END DO 459 ! III - Barotropic repartition of the sea surface height over the baroclinic profile 460 ! ================================================================================== 461 ! add ( ssh increment + "baroclinicity error" ) proportionly to e3t(n) 462 ! - ML - baroclinicity error should be better treated in the future 463 ! i.e. locally and not spread over the water column. 464 ! (keep in mind that the idea is to reduce Eulerian velocity as much as possible) 465 zht(:,:) = 0. 627 zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * (hdivn(:,:,jk) - hdivb(:,:,jk)) 628 END DO 629 zhdiv(:,:) = zhdiv(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) ) 630 631 IF( neuler == 0 .AND. kt == nit000 ) THEN 632 z2dt = rdt 633 ELSE 634 z2dt = 2.0_wp * rdt 635 ENDIF 636 466 637 DO jk = 1, jpkm1 467 zht(:,:) = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk) 468 END DO 469 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 638 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - z2dt * fse3t_n(:,:,jk) * & 639 & (hdivn(:,:,jk) - hdivb(:,:,jk) - zhdiv(:,:)) 640 END DO 641 ENDIF 642 643 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde or layer coordinate ! 644 ! ! ---baroclinic part--------- ! 470 645 DO jk = 1, jpkm1 471 dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 472 END DO 473 646 fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + (tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk)) 647 END DO 474 648 ENDIF 475 649 476 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde or layer coordinate ! 477 ! ! ---baroclinic part--------- ! 650 IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN ! - ML - test: control prints for debuging 651 ! 652 zht(:,:) = 0.0_wp 478 653 DO jk = 1, jpkm1 479 fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk) 480 END DO 481 ENDIF 482 483 IF( ln_vvl_dbg .AND. .NOT. ll_do_bclinic ) THEN ! - ML - test: control prints for debuging 484 ! 654 zht(:,:) = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk) 655 END DO 485 656 IF( lwp ) WRITE(numout, *) 'kt =', kt 486 657 IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 487 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) )658 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ), mask = tmask(:,:,1) == 1.e0 ) 488 659 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 489 660 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax 490 661 END IF 491 662 ! 663 z_tmin = MINVAL( fse3t_n(:,:,:), mask = tmask(:,:,:) == 1.e0 ) 664 IF( lk_mpp ) CALL mpp_min( z_tmin ) ! min over the global domain 665 IF( lwp ) WRITE(numout, *) kt,' MINVAL(fse3t_n) =', z_tmin 666 ! 667 z_tmin = MINVAL( fse3u_n(:,:,:), mask = umask(:,:,:) == 1.e0 ) 668 IF( lk_mpp ) CALL mpp_min( z_tmin ) ! min over the global domain 669 IF( lwp ) WRITE(numout, *) kt,' MINVAL(fse3u_n) =', z_tmin 670 ! 671 z_tmin = MINVAL( fse3v_n(:,:,:), mask = vmask(:,:,:) == 1.e0 ) 672 IF( lk_mpp ) CALL mpp_min( z_tmin ) ! min over the global domain 673 IF( lwp ) WRITE(numout, *) kt,' MINVAL(fse3v_n) =', z_tmin 674 ! 675 z_tmin = MINVAL( fse3f_n(:,:,:), mask = fmask(:,:,:) == 1.e0 ) 676 IF( lk_mpp ) CALL mpp_min( z_tmin ) ! min over the global domain 677 IF( lwp ) WRITE(numout, *) kt,' MINVAL(fse3f_n) =', z_tmin 678 ! 492 679 zht(:,:) = 0.0_wp 493 680 DO jk = 1, jpkm1 … … 496 683 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) ) 497 684 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 498 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn -SUM(fse3t_n))) =', z_tmax685 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn -SUM(fse3t_n))) =', z_tmax 499 686 ! 500 687 zht(:,:) = 0.0_wp 501 688 DO jk = 1, jpkm1 502 zht(:,:) = zht(:,:) + fse3t_a(:,:,jk) * tmask(:,:,jk) 503 END DO 504 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) ) 689 zht(:,:) = zht(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk) 690 END DO 691 zwu(:,:) = 0._wp 692 DO jj = 1, jpjm1 693 DO ji = 1, fs_jpim1 ! vector opt. 694 zwu(ji,jj) = 0.5_wp * umask(ji,jj,1) * r1_e12u(ji,jj) & 695 & * ( e12t(ji,jj) * sshn(ji,jj) + e12t(ji+1,jj) * sshn(ji+1,jj) ) 696 END DO 697 END DO 698 CALL lbc_lnk( zwu(:,:), 'U', 1._wp ) 699 z_tmax = MAXVAL( umask(:,:,1) * umask_i(:,:) * ABS( hu_0(:,:) + zwu(:,:) - zht(:,:) ) ) 505 700 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 506 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(h t_0+ssha-SUM(fse3t_a))) =', z_tmax701 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(hu_0+sshu_n-SUM(fse3u_n))) =', z_tmax 507 702 ! 508 703 zht(:,:) = 0.0_wp 509 704 DO jk = 1, jpkm1 510 zht(:,:) = zht(:,:) + fse3t_b(:,:,jk) * tmask(:,:,jk) 511 END DO 512 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) ) 705 zht(:,:) = zht(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk) 706 END DO 707 zwv(:,:) = 0._wp 708 DO jj = 1, jpjm1 709 DO ji = 1, fs_jpim1 ! vector opt. 710 zwv(ji,jj) = 0.5_wp * vmask(ji,jj,1) * r1_e12v(ji,jj) & 711 & * ( e12t(ji,jj) * sshn(ji,jj) + e12t(ji,jj+1) * sshn(ji,jj+1) ) 712 END DO 713 END DO 714 CALL lbc_lnk( zwv(:,:), 'V', 1._wp ) 715 z_tmax = MAXVAL( vmask(:,:,1) * vmask_i(:,:) * ABS( hv_0(:,:) + zwv(:,:) - zht(:,:) ) ) 513 716 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 514 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(fse3t_b))) =', z_tmax 515 ! 516 z_tmax = MAXVAL( tmask(:,:,1) * ABS( sshb(:,:) ) ) 717 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(hv_0+sshv_n-SUM(fse3v_n))) =', z_tmax 718 ! 719 zht(:,:) = 0.0_wp 720 DO jk = 1, jpkm1 721 DO jj = 1, jpjm1 722 zht(:,jj) = zht(:,jj) + fse3f_n(:,jj,jk) * umask(:,jj,jk)*umask(:,jj+1,jk) 723 END DO 724 END DO 725 zwu(:,:) = 0._wp 726 DO jj = 1, jpjm1 727 DO ji = 1, fs_jpim1 ! vector opt. 728 zwu(ji,jj) = 0.25_wp * umask(ji,jj,1) * umask(ji,jj+1,1) * r1_e12f(ji,jj) & 729 & * ( e12t(ji ,jj) * sshn(ji ,jj) + e12t(ji ,jj+1) * sshn(ji ,jj+1) & 730 & + e12t(ji+1,jj) * sshn(ji+1,jj) + e12t(ji+1,jj+1) * sshn(ji+1,jj+1) ) 731 END DO 732 END DO 733 CALL lbc_lnk( zht(:,:), 'F', 1._wp ) 734 CALL lbc_lnk( zwu(:,:), 'F', 1._wp ) 735 z_tmax = MAXVAL( fmask(:,:,1) * fmask_i(:,:) * ABS( hf_0(:,:) + zwu(:,:) - zht(:,:) ) ) 517 736 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 518 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(sshb))) =', z_tmax 519 ! 520 z_tmax = MAXVAL( tmask(:,:,1) * ABS( sshn(:,:) ) ) 521 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 522 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(sshn))) =', z_tmax 523 ! 524 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssha(:,:) ) ) 525 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 526 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssha))) =', z_tmax 737 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(hf_0+sshf_n-SUM(fse3f_n))) =', z_tmax 738 ! 527 739 END IF 528 740 … … 549 761 550 762 CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv ) 551 CALL wrk_dealloc( jpi, jpj, jpk, ze3t 763 CALL wrk_dealloc( jpi, jpj, jpk, ze3t, ztu, ztv ) 552 764 553 765 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_sf_nxt') … … 583 795 INTEGER, INTENT( in ) :: kt ! time step 584 796 !! * Local declarations 585 INTEGER :: ji,jj,jk ! dummy loop indices586 REAL(wp) :: zcoef797 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwrk 798 INTEGER :: jk ! dummy loop indices 587 799 !!---------------------------------------------------------------------- 588 800 … … 607 819 tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 608 820 ENDIF 821 609 822 fsdept_b(:,:,:) = fsdept_n(:,:,:) 610 823 fsdepw_b(:,:,:) = fsdepw_n(:,:,:) … … 620 833 ! - ML - fse3u_b and fse3v_b are allready computed in dynnxt 621 834 ! - JC - hu_b, hv_b, hur_b, hvr_b also 622 CALL dom_vvl_interpol( fse3 u_n(:,:,:), fse3f_n (:,:,:), 'F' )835 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3f_n (:,:,:), 'F' ) 623 836 ! Vertical scale factor interpolations 624 837 ! ------------------------------------ … … 631 844 ! t- and w- points depth 632 845 ! ---------------------- 633 ! set the isf depth as it is in the initial step634 846 fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 635 847 fsdepw_n(:,:,1) = 0.0_wp 636 848 fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 637 638 849 DO jk = 2, jpk 639 DO jj = 1,jpj 640 DO ji = 1,jpi 641 ! zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) ! 0 everywhere tmask = wmask, ie everywhere expect at jk = mikt 642 ! 1 for jk = mikt 643 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 644 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 645 fsdept_n(ji,jj,jk) = zcoef * ( fsdepw_n(ji,jj,jk ) + 0.5 * fse3w_n(ji,jj,jk)) & 646 & + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)) 647 fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - sshn(ji,jj) 648 END DO 649 END DO 850 fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk) 851 fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1) 852 fsde3w_n(:,:,jk) = fsdept_n(:,:,jk ) - sshn (:,:) 650 853 END DO 651 652 854 ! Local depth and Inverse of the local depth of the water column at u- and v- points 653 855 ! ---------------------------------------------------------------------------------- … … 666 868 END DO 667 869 870 ! Write additional diagnostics 871 ! ============================ 872 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) CALL dom_vvl_dia( kt) 873 668 874 ! write restart file 669 875 ! ================== … … 674 880 END SUBROUTINE dom_vvl_sf_swp 675 881 882 SUBROUTINE dom_vvl_dia( kt ) 883 !!---------------------------------------------------------------------- 884 !! *** ROUTINE dom_vvl_dia *** 885 !! 886 !! ** Purpose : Output some diagnostics in ztilde/zlayer cases 887 !! 888 !!---------------------------------------------------------------------- 889 !! * Arguments 890 INTEGER, INTENT( in ) :: kt ! time step 891 !! * Local declarations 892 INTEGER :: ji,jj,jk ! dummy loop indices 893 REAL(wp) :: zufim1, zufi, zvfjm1, zvfj, ztmp1, z2dt 894 REAL(wp), DIMENSION(4) :: zr1 895 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwdw, zout 896 !!---------------------------------------------------------------------- 897 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_dia') 898 ! 899 CALL wrk_alloc( jpi, jpj, jpk, zwdw, zout ) 900 ! 901 ! Compute internal interfaces depths: 902 !------------------------------------ 903 IF ( iom_use("dh_tilde").OR.iom_use("depw_tilde").OR.iom_use("stiff_tilde")) THEN 904 zwdw(:,:,1) = 0.e0 905 DO jj = 1, jpj 906 DO ji = 1, jpi 907 DO jk = 2, jpkm1 908 zwdw(ji,jj,jk) = zwdw(ji,jj,jk-1) + & 909 & (tilde_e3t_n(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1) 910 END DO 911 END DO 912 END DO 913 ENDIF 914 ! 915 ! Output interface depth anomaly: 916 ! ------------------------------- 917 IF ( iom_use("depw_tilde") ) CALL iom_put( "depw_tilde", (zwdw(:,:,:)-gdepw_0(:,:,:))*tmask(:,:,:) ) 918 ! 919 ! Output grid stiffness (T-points): 920 ! --------------------------------- 921 IF ( iom_use("stiff_tilde" ) ) THEN 922 zr1(:) = 0.e0 923 zout(:,:,jpk) = 0.e0 924 DO ji = 2, jpim1 925 DO jj = 2, jpjm1 926 DO jk = 1, jpkm1 927 zr1(1) = umask(ji-1,jj ,jk) *abs( (zwdw(ji ,jj ,jk )-zwdw(ji-1,jj ,jk ) & 928 & +zwdw(ji ,jj ,jk+1)-zwdw(ji-1,jj ,jk+1)) & 929 & /(zwdw(ji ,jj ,jk )+zwdw(ji-1,jj ,jk ) & 930 & -zwdw(ji ,jj ,jk+1)-zwdw(ji-1,jj ,jk+1) + rsmall) ) 931 zr1(2) = umask(ji ,jj ,jk) *abs( (zwdw(ji+1,jj ,jk )-zwdw(ji ,jj ,jk ) & 932 & +zwdw(ji+1,jj ,jk+1)-zwdw(ji ,jj ,jk+1)) & 933 & /(zwdw(ji+1,jj ,jk )+zwdw(ji ,jj ,jk ) & 934 & -zwdw(ji+1,jj ,jk+1)-zwdw(ji ,jj ,jk+1) + rsmall) ) 935 zr1(3) = vmask(ji ,jj ,jk) *abs( (zwdw(ji ,jj+1,jk )-zwdw(ji ,jj ,jk ) & 936 & +zwdw(ji ,jj+1,jk+1)-zwdw(ji ,jj ,jk+1)) & 937 & /(zwdw(ji ,jj+1,jk )+zwdw(ji ,jj ,jk ) & 938 & -zwdw(ji ,jj+1,jk+1)-zwdw(ji ,jj ,jk+1) + rsmall) ) 939 zr1(4) = vmask(ji ,jj-1,jk) *abs( (zwdw(ji ,jj ,jk )-zwdw(ji ,jj-1,jk ) & 940 & +zwdw(ji ,jj ,jk+1)-zwdw(ji ,jj-1,jk+1)) & 941 & /(zwdw(ji ,jj ,jk )+zwdw(ji ,jj-1,jk ) & 942 & -zwdw(ji, jj ,jk+1)-zwdw(ji ,jj-1,jk+1) + rsmall) ) 943 zout(ji,jj,jk) = MAXVAL(zr1(1:4)) 944 END DO 945 END DO 946 END DO 947 948 CALL lbc_lnk( zout, 'T', 1. ) 949 CALL iom_put( "stiff_tilde", zout(:,:,:) ) 950 END IF 951 ! Output Horizontal Laplacian of interfaces depths (W-points): 952 ! ------------------------------------------------------------ 953 IF ( iom_use("dh_tilde") ) THEN 954 ! 955 zout(:,:,1 )=0._wp 956 zout(:,:,jpk)=0._wp 957 DO jk = 2, jpkm1 958 DO jj = 2, jpjm1 959 DO ji = fs_2, fs_jpim1 ! vector opt. 960 zufim1 = umask(ji-1,jj,jk) * re2u_e1u(ji-1,jj) * ( zwdw(ji-1,jj,jk) - zwdw(ji ,jj ,jk) ) 961 zufi = umask(ji ,jj,jk) * re2u_e1u(ji ,jj) * ( zwdw(ji ,jj,jk) - zwdw(ji+1,jj ,jk) ) 962 zvfjm1 = vmask(ji,jj-1,jk) * re1v_e2v(ji,jj-1) * ( zwdw(ji,jj-1,jk) - zwdw(ji ,jj ,jk) ) 963 zvfj = vmask(ji,jj ,jk) * re1v_e2v(ji,jj ) * ( zwdw(ji,jj ,jk) - zwdw(ji ,jj+1,jk) ) 964 ! ztmp1 = (zufim1-zufi+zvfjm1-zvfj) * SQRT(r1_e12t(ji,jj)) 965 ztmp1 = (zufim1-zufi+zvfjm1-zvfj) * r1_e12t(ji,jj) 966 ! zout(ji,jj,jk) = ABS(ztmp1)*tmask(ji,jj,jk) 967 zout(ji,jj,jk) = ztmp1 968 END DO 969 END DO 970 END DO 971 ! Mask open boundaries: 972 #if defined key_bdy 973 IF (lk_bdy) THEN 974 DO jk = 1, jpkm1 975 zout(:,:,jk) = zout(:,:,jk) * bdytmask(:,:) 976 END DO 977 ENDIF 978 #endif 979 CALL lbc_lnk( zout(:,:,:), 'T', 1. ) 980 zwdw(:,:,:) = zout(:,:,:) 981 DO jk = 2, jpkm1 982 DO jj = 2, jpjm1 983 DO ji = fs_2, fs_jpim1 ! vector opt. 984 zufim1 = umask(ji-1,jj,jk) * re2u_e1u(ji-1,jj) * ( zwdw(ji-1,jj,jk) - zwdw(ji ,jj ,jk) ) 985 zufi = umask(ji ,jj,jk) * re2u_e1u(ji ,jj) * ( zwdw(ji ,jj,jk) - zwdw(ji+1,jj ,jk) ) 986 zvfjm1 = vmask(ji,jj-1,jk) * re1v_e2v(ji,jj-1) * ( zwdw(ji,jj-1,jk) - zwdw(ji ,jj ,jk) ) 987 zvfj = vmask(ji,jj ,jk) * re1v_e2v(ji,jj ) * ( zwdw(ji,jj ,jk) - zwdw(ji ,jj+1,jk) ) 988 ! ztmp1 = (zufim1-zufi+zvfjm1-zvfj) * SQRT(r1_e12t(ji,jj)) 989 ztmp1 = (zufim1-zufi+zvfjm1-zvfj) * r1_e12t(ji,jj) 990 zout(ji,jj,jk) = ABS(ztmp1)*tmask(ji,jj,jk) 991 END DO 992 END DO 993 END DO 994 ! Mask open boundaries: 995 #if defined key_bdy 996 IF (lk_bdy) THEN 997 DO jk = 1, jpkm1 998 zout(:,:,jk) = zout(:,:,jk) * bdytmask(:,:) 999 END DO 1000 ENDIF 1001 #endif 1002 CALL lbc_lnk( zout(:,:,:), 'T', 1. ) 1003 ! 1004 CALL iom_put( "dh_tilde", zout(:,:,:) ) 1005 ENDIF 1006 ! 1007 ! Output vertical Laplacian of interfaces depths (W-points): 1008 ! ---------------------------------------------------------- 1009 IF ( iom_use("dz_tilde" ) ) THEN 1010 zout(:,:,1 ) = 0._wp 1011 zout(:,:,jpk) = 0._wp 1012 DO jk=2,jpkm1 1013 zout(:,:,jk) = 2._wp*ABS(tilde_e3t_n(:,:,jk)+e3t_0(:,:,jk)-tilde_e3t_n(:,:,jk-1)-e3t_0(:,:,jk-1)) & 1014 & /(tilde_e3t_n(:,:,jk)+e3t_0(:,:,jk)+tilde_e3t_n(:,:,jk-1)+e3t_0(:,:,jk-1)) & 1015 & * tmask(:,:,jk) 1016 END DO 1017 CALL iom_put( "dz_tilde", zout(:,:,:) ) 1018 END IF 1019 ! 1020 ! 1021 ! Output low pass U-velocity: 1022 ! --------------------------- 1023 IF ( iom_use("un_lf_tilde" ).AND.ln_vvl_ztilde ) THEN 1024 zout(:,:,jpk) = 0.e0 1025 DO jk=1,jpkm1 1026 zout(:,:,jk) = un_lf(:,:,jk)/fse3u_n(:,:,jk)*r1_e2u(:,:) 1027 END DO 1028 CALL iom_put( "un_lf_tilde", zout(:,:,:) ) 1029 END IF 1030 ! 1031 ! Output low pass V-velocity: 1032 ! --------------------------- 1033 IF ( iom_use("vn_lf_tilde" ).AND.ln_vvl_ztilde ) THEN 1034 zout(:,:,jpk) = 0.e0 1035 DO jk=1,jpkm1 1036 zout(:,:,jk) = vn_lf(:,:,jk)/fse3v_n(:,:,jk)*r1_e1v(:,:) 1037 END DO 1038 CALL iom_put( "vn_lf_tilde", zout(:,:,:) ) 1039 END IF 1040 ! 1041 CALL wrk_dealloc( jpi, jpj, jpk, zwdw, zout ) 1042 ! 1043 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_dia') 1044 ! 1045 END SUBROUTINE dom_vvl_dia 676 1046 677 1047 SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout ) … … 691 1061 ! ! = 'U', 'V', 'W, 'F', 'UW' or 'VW' 692 1062 !! * Local declarations 1063 INTEGER :: nmet ! horizontal interpolation method 693 1064 INTEGER :: ji, jj, jk ! dummy loop indices 1065 INTEGER :: jkbot ! bottom level 694 1066 LOGICAL :: l_is_orca ! local logical 1067 REAL(wp) :: zmin, zdo, zup, ztap 1068 REAL(wp), POINTER, DIMENSION(:,:) :: zs ! Surface interface depth anomaly 1069 REAL(wp), POINTER, DIMENSION(:,:,:) :: zw ! Interface depth anomaly 695 1070 !!---------------------------------------------------------------------- 696 1071 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_interpol') … … 698 1073 l_is_orca = .FALSE. 699 1074 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) l_is_orca = .TRUE. ! ORCA R2 configuration - will need to correct some locations 1075 1076 1077 nmet=1 ! Original method (Surely wrong) 1078 ! nmet=1 ! Interface interpolation 1079 ! nmet=2 ! Internal interfaces interpolation only, spread barotropic increment 1080 ! nmet=1 1081 ztap = 0.1_wp ! Minimum fraction of T-point thickness at cell interfaces 1082 1083 IF ( (nmet==1).OR.(nmet==2) ) THEN 1084 SELECT CASE ( pout ) 1085 ! 1086 CASE( 'U', 'V', 'F' ) 1087 ! Compute interface depth anomaly at T-points 1088 CALL wrk_alloc( jpi, jpj, jpk, zw ) 1089 CALL wrk_alloc( jpi, jpj, zs ) 1090 ! 1091 zw(:,:,:) = 0._wp 1092 ! 1093 DO jj = 1, jpj 1094 DO ji = 1, jpi 1095 jkbot = mbkt(ji,jj) 1096 DO jk=jkbot,1,-1 1097 zw(ji,jj,jk) = zw(ji,jj,jk+1) - pe3_in(ji,jj,jk) + e3t_0(ji,jj,jk) 1098 END DO 1099 END DO 1100 END DO 1101 ! 1102 ! 1103 IF (nmet==2) THEN ! Consider "internal" interfaces only 1104 zs(:,:) = - zw(:,:,1) ! Save surface anomaly (ssh) 1105 ! 1106 zw(:,:,:) = 0._wp 1107 DO jj = 1, jpj 1108 DO ji = 1, jpi 1109 jkbot = mbkt(ji,jj) 1110 DO jk=jkbot,1,-1 1111 zw(ji,jj,jk) = zw(ji,jj,jk+1) - ( pe3_in(ji,jj,jk) & 1112 & * ht_0(ji,jj) / (ht_0(ji,jj) + zs(ji,jj) + 1._wp - tmask(ji,jj,1)) & 1113 & - e3t_0(ji,jj,jk)) * tmask(ji,jj,jk) 1114 END DO 1115 END DO 1116 END DO 1117 ENDIF 1118 ! 1119 END SELECT 1120 END IF 1121 1122 pe3_out(:,:,:) = 0._wp 700 1123 701 1124 SELECT CASE ( pout ) 702 1125 ! ! ------------------------------------- ! 703 1126 CASE( 'U' ) ! interpolation from T-point to U-point ! 704 ! ! ------------------------------------- ! 1127 ! ! ------------------------------------- ! 705 1128 ! horizontal surface weighted interpolation 706 DO jk = 1, jpk 1129 IF (nmet==0) THEN 1130 DO jk = 1, jpk 1131 DO jj = 1, jpjm1 1132 DO ji = 1, fs_jpim1 ! vector opt. 1133 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj) & 1134 & * ( e12t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & 1135 & + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 1136 END DO 1137 END DO 1138 END DO 1139 ENDIF 1140 ! 1141 IF ( (nmet==1).OR.(nmet==2) ) THEN 707 1142 DO jj = 1, jpjm1 708 1143 DO ji = 1, fs_jpim1 ! vector opt. 709 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj) & 710 & * ( e12t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & 711 & + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 712 END DO 713 END DO 714 END DO 1144 ! Correction at last level: 1145 jkbot = mbku(ji,jj) 1146 pe3_out(ji,jj,jkbot) = -0.5_wp * umask(ji,jj,jkbot) * r1_e12u(ji,jj) & 1147 & * ( e12t(ji ,jj) * zw(ji ,jj,jkbot) & 1148 & + e12t(ji+1,jj) * zw(ji+1,jj,jkbot) ) 1149 ! 1150 ! If there is a step, taper bottom interface: 1151 IF (hu_0(ji,jj) < 0.5_wp * ( ht_0(ji,jj) + ht_0(ji+1,jj) ) ) THEN 1152 IF ( ht_0(ji+1,jj) < ht_0(ji,jj) ) THEN 1153 ! zmin = ztap * pe3_in(ji+1,jj,jkbot) 1154 zmin = ztap * (-zw(ji+1,jj,jkbot)+e3t_0(ji+1,jj,jkbot)) 1155 ELSE 1156 ! zmin = ztap * pe3_in(ji ,jj,jkbot) 1157 zmin = ztap * (-zw(ji,jj,jkbot)+e3t_0(ji,jj,jkbot)) 1158 ENDIF 1159 zmin = -e3u_0(ji,jj,jkbot) + zmin 1160 pe3_out(ji,jj,jkbot) = MAX(pe3_out(ji,jj,jkbot), zmin) 1161 ENDIF 1162 ! 1163 zdo = -pe3_out(ji,jj,jkbot) 1164 DO jk=jkbot-1,1,-1 1165 zup = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji ,jj) & 1166 & *( e12t(ji ,jj) * zw(ji ,jj,jk) & 1167 & +e12t(ji+1,jj) * zw(ji+1,jj,jk) ) 1168 pe3_out(ji,jj,jk) = zdo - zup 1169 zdo = zdo - pe3_out(ji,jj,jk) 1170 END DO 1171 END DO 1172 END DO 1173 END IF 1174 ! 1175 IF (nmet==2) THEN ! Spread sea level anomaly 1176 DO jj = 1, jpjm1 1177 DO ji = 1, fs_jpim1 ! vector opt. 1178 DO jk=1,jpk 1179 pe3_out(ji,jj,jk) = pe3_out(ji,jj,jk) & 1180 & + ( pe3_out(ji,jj,jk) + e3u_0(ji,jj,jk) ) & 1181 & / ( hu_0(ji,jj) + 1._wp - umask_i(ji,jj) ) & 1182 & * 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj) & 1183 & * ( e12t(ji,jj)*zs(ji,jj) + e12t(ji+1,jj)*zs(ji+1,jj) ) 1184 END DO 1185 END DO 1186 END DO 1187 ! 1188 ENDIF 715 1189 ! 716 1190 IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout ) … … 718 1192 CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp ) 719 1193 pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:) 1194 ! 1195 IF ( (nmet==1).OR.(nmet==2) ) CALL wrk_dealloc( jpi, jpj, zs ) 1196 IF ( (nmet==1).OR.(nmet==2) ) CALL wrk_dealloc( jpi, jpj, jpk, zw ) 1197 ! 720 1198 ! ! ------------------------------------- ! 721 1199 CASE( 'V' ) ! interpolation from T-point to V-point ! 722 1200 ! ! ------------------------------------- ! 723 1201 ! horizontal surface weighted interpolation 724 DO jk = 1, jpk 1202 IF (nmet==0) THEN 1203 DO jk = 1, jpk 1204 DO jj = 1, jpjm1 1205 DO ji = 1, fs_jpim1 ! vector opt. 1206 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj) & 1207 & * ( e12t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & 1208 & + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 1209 END DO 1210 END DO 1211 END DO 1212 ENDIF 1213 ! 1214 IF ( (nmet==1).OR.(nmet==2) ) THEN 725 1215 DO jj = 1, jpjm1 726 1216 DO ji = 1, fs_jpim1 ! vector opt. 727 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj) & 728 & * ( e12t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & 729 & + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 730 END DO 731 END DO 732 END DO 1217 ! Correction at last level: 1218 jkbot = mbkv(ji,jj) 1219 pe3_out(ji,jj,jkbot) = -0.5_wp * vmask(ji,jj,jkbot) * r1_e12v(ji,jj) & 1220 & * ( e12t(ji,jj ) * zw(ji,jj ,jkbot) & 1221 & + e12t(ji,jj+1) * zw(ji,jj+1,jkbot) ) 1222 ! 1223 ! If there is a step, taper bottom interface: 1224 IF (hv_0(ji,jj) < 0.5_wp * ( ht_0(ji,jj) + ht_0(ji,jj+1) ) ) THEN 1225 IF ( ht_0(ji,jj+1) < ht_0(ji,jj) ) THEN 1226 ! zmin = ztap * pe3_in(ji,jj+1,jkbot) 1227 zmin = ztap * (-zw(ji,jj+1,jkbot)+e3t_0(ji,jj+1,jkbot)) 1228 ELSE 1229 ! zmin = ztap * pe3_in(ji,jj ,jkbot) 1230 zmin = ztap * (-zw(ji,jj,jkbot)+e3t_0(ji,jj,jkbot)) 1231 ENDIF 1232 zmin = -e3v_0(ji,jj,jkbot) + zmin 1233 pe3_out(ji,jj,jkbot) = MAX(pe3_out(ji,jj,jkbot), zmin) 1234 ENDIF 1235 ! 1236 zdo = -pe3_out(ji,jj,jkbot) 1237 DO jk=jkbot-1,1,-1 1238 zup = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj ) & 1239 & * ( e12t(ji,jj ) * zw(ji,jj ,jk) & 1240 & +e12t(ji,jj+1) * zw(ji,jj+1,jk) ) 1241 ! 1242 pe3_out(ji,jj,jk) = zdo - zup 1243 zdo = zdo - pe3_out(ji,jj,jk) 1244 END DO 1245 END DO 1246 END DO 1247 END IF 1248 ! 1249 IF (nmet==2) THEN ! Spread sea level anomaly 1250 ! 1251 DO jj = 1, jpjm1 1252 DO ji = 1, fs_jpim1 ! vector opt. 1253 DO jk=1,jpk 1254 pe3_out(ji,jj,jk) = pe3_out(ji,jj,jk) & 1255 & + ( pe3_out(ji,jj,jk) + e3v_0(ji,jj,jk) ) & 1256 & / ( hv_0(ji,jj) + 1._wp - vmask_i(ji,jj) ) & 1257 & * 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj) & 1258 & * ( e12t(ji,jj)*zs(ji,jj) + e12t(ji,jj+1)*zs(ji,jj+1) ) 1259 END DO 1260 END DO 1261 END DO 1262 ! 1263 ENDIF 733 1264 ! 734 1265 IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout ) … … 736 1267 CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp ) 737 1268 pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:) 1269 ! 1270 IF ( (nmet==1).OR.(nmet==2) ) CALL wrk_dealloc( jpi, jpj, zs ) 1271 IF ( (nmet==1).OR.(nmet==2) ) CALL wrk_dealloc( jpi, jpj, jpk, zw ) 1272 ! 738 1273 ! ! ------------------------------------- ! 739 CASE( 'F' ) ! interpolation from U-point to F-point !1274 CASE( 'F' ) ! interpolation from T-point to F-point ! 740 1275 ! ! ------------------------------------- ! 741 1276 ! horizontal surface weighted interpolation 742 DO jk = 1, jpk 1277 IF (nmet==0) THEN 1278 DO jk=1,jpk 1279 DO jj = 1, jpjm1 1280 DO ji = 1, fs_jpim1 ! vector opt. 1281 pe3_out(ji,jj,jk) = 0.25_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj) & 1282 & * ( e12t(ji ,jj ) * ( pe3_in(ji ,jj ,jk)-e3t_0(ji ,jj ,jk) ) & 1283 & + e12t(ji ,jj+1) * ( pe3_in(ji ,jj+1,jk)-e3t_0(ji ,jj+1,jk) ) & 1284 & + e12t(ji+1,jj ) * ( pe3_in(ji+1,jj ,jk)-e3t_0(ji+1,jj ,jk) ) & 1285 & + e12t(ji+1,jj+1) * ( pe3_in(ji+1,jj+1,jk)-e3t_0(ji+1,jj+1,jk) )) 1286 END DO 1287 END DO 1288 END DO 1289 ENDIF 1290 ! 1291 IF ( (nmet==1).OR.(nmet==2) ) THEN 743 1292 DO jj = 1, jpjm1 744 1293 DO ji = 1, fs_jpim1 ! vector opt. 745 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj) & 746 & * ( e12u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) & 747 & + e12u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 748 END DO 749 END DO 750 END DO 1294 ! bottom correction: 1295 jkbot = MIN(mbku(ji,jj), mbku(ji,jj+1)) 1296 pe3_out(ji,jj,jkbot) = -0.25_wp * umask(ji,jj,jkbot) * umask(ji,jj+1,jkbot) * r1_e12f(ji,jj) & 1297 & * ( e12t(ji ,jj ) * zw(ji ,jj ,jkbot) & 1298 & + e12t(ji+1,jj ) * zw(ji+1,jj ,jkbot) & 1299 & + e12t(ji ,jj+1) * zw(ji ,jj+1,jkbot) & 1300 & + e12t(ji+1,jj+1) * zw(ji+1,jj+1,jkbot) ) 1301 ! 1302 ! If there is a step, taper bottom interface: 1303 IF (hf_0(ji,jj) < 0.5_wp * ( hu_0(ji,jj ) + hu_0(ji,jj+1) ) ) THEN 1304 IF ( hu_0(ji,jj+1) < hu_0(ji,jj) ) THEN 1305 IF ( ht_0(ji+1,jj+1) < ht_0(ji ,jj+1) ) THEN 1306 ! zmin = ztap * pe3_in(ji+1,jj+1,jkbot) 1307 zmin = ztap * (-zw(ji+1,jj+1,jkbot)+e3t_0(ji+1,jj+1,jkbot)) 1308 ELSE 1309 ! zmin = ztap * pe3_in(ji ,jj+1,jkbot) 1310 zmin = ztap * (-zw(ji,jj+1,jkbot)+e3t_0(ji,jj+1,jkbot)) 1311 ENDIF 1312 ELSE 1313 IF ( ht_0(ji+1,jj ) < ht_0(ji ,jj ) ) THEN 1314 ! zmin = ztap * pe3_in(ji+1,jj ,jkbot) 1315 zmin = ztap * (-zw(ji+1,jj,jkbot)+e3t_0(ji+1,jj,jkbot)) 1316 ELSE 1317 ! zmin = ztap * pe3_in(ji ,jj ,jkbot) 1318 zmin = ztap * (-zw(ji,jj,jkbot)+e3t_0(ji,jj,jkbot)) 1319 ENDIF 1320 ENDIF 1321 zmin = -e3f_0(ji,jj,jkbot) + zmin 1322 pe3_out(ji,jj,jkbot) = MAX(pe3_out(ji,jj,jkbot), zmin) 1323 ENDIF 1324 ! 1325 zdo = -pe3_out(ji,jj,jkbot) 1326 DO jk=jkbot-1,1,-1 1327 zup = 0.25_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj) & 1328 & * ( e12t(ji ,jj ) * zw(ji ,jj ,jk) & 1329 & + e12t(ji+1,jj ) * zw(ji+1,jj ,jk) & 1330 & + e12t(ji ,jj+1) * zw(ji ,jj+1,jk) & 1331 & + e12t(ji+1,jj+1) * zw(ji+1,jj+1,jk) ) 1332 pe3_out(ji,jj,jk) = zdo - zup 1333 ! 1334 zdo = zdo - pe3_out(ji,jj,jk) 1335 END DO 1336 ! 1337 END DO 1338 END DO 1339 ENDIF 1340 ! 1341 IF (nmet==2) THEN ! Spread sea level anomaly 1342 ! 1343 DO jj = 1, jpjm1 1344 DO ji = 1, fs_jpim1 ! vector opt. 1345 DO jk=1,jpk 1346 pe3_out(ji,jj,jk) = pe3_out(ji,jj,jk) & 1347 & + ( pe3_out(ji,jj,jk) + e3f_0(ji,jj,jk) ) & 1348 & / ( hf_0(ji,jj) + 1._wp - umask_i(ji,jj)*umask_i(ji,jj+1) ) & 1349 & * 0.25_wp * umask(ji,jj,jk)*umask(ji,jj+1,jk)*r1_e12f(ji,jj) & 1350 & * ( e12t(ji ,jj)*zs(ji ,jj) + e12t(ji ,jj+1)*zs(ji ,jj+1) & 1351 & +e12t(ji+1,jj)*zs(ji+1,jj) + e12t(ji+1,jj+1)*zs(ji+1,jj+1) ) 1352 END DO 1353 END DO 1354 END DO 1355 END IF 751 1356 ! 752 1357 IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout ) … … 754 1359 CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp ) 755 1360 pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:) 1361 ! 1362 IF ( (nmet==1).OR.(nmet==2) ) CALL wrk_dealloc( jpi, jpj, zs ) 1363 IF ( (nmet==1).OR.(nmet==2) ) CALL wrk_dealloc( jpi, jpj, jpk, zw ) 1364 ! 756 1365 ! ! ------------------------------------- ! 757 1366 CASE( 'W' ) ! interpolation from T-point to W-point ! … … 808 1417 !! * Local declarations 809 1418 INTEGER :: jk 810 INTEGER :: id1, id2, id3, id4, id5 ! local integers 1419 INTEGER, DIMENSION(7) :: id ! local integers 1420 REAL(wp), POINTER, DIMENSION(:,:) :: zhdiv 811 1421 !!---------------------------------------------------------------------- 812 1422 ! … … 818 1428 CALL iom_get( numror, jpdom_autoglo, 'sshn' , sshn ) 819 1429 ! 820 id1 = iom_varid( numror, 'fse3t_b', ldstop = .FALSE. ) 821 id2 = iom_varid( numror, 'fse3t_n', ldstop = .FALSE. ) 822 id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. ) 823 id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) 824 id5 = iom_varid( numror, 'hdiv_lf', ldstop = .FALSE. ) 1430 id(1) = iom_varid( numror, 'fse3t_b', ldstop = .FALSE. ) 1431 id(2) = iom_varid( numror, 'fse3t_n', ldstop = .FALSE. ) 1432 id(3) = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. ) 1433 id(4) = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) 1434 id(5) = iom_varid( numror, 'hdivn_lf', ldstop = .FALSE. ) 1435 id(6) = iom_varid( numror, 'un_lf' , ldstop = .FALSE. ) 1436 id(7) = iom_varid( numror, 'vn_lf' , ldstop = .FALSE. ) 1437 825 1438 ! ! --------- ! 826 1439 ! ! all cases ! 827 1440 ! ! --------- ! 828 IF( MIN( id 1, id2) > 0 ) THEN ! all required arrays exist1441 IF( MIN( id(1), id(2) ) > 0 ) THEN ! all required arrays exist 829 1442 CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) ) 830 1443 CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) ) … … 838 1451 fse3t_b(:,:,:) = fse3t_n(:,:,:) 839 1452 ENDIF 840 ELSE IF( id 1> 0 ) THEN1453 ELSE IF( id(1) > 0 ) THEN 841 1454 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart files' 842 1455 IF(lwp) write(numout,*) 'fse3t_n set equal to fse3t_b.' 843 1456 IF(lwp) write(numout,*) 'neuler is forced to 0' 844 1457 CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) ) 845 fse3t_ n(:,:,:) = fse3t_b(:,:,:)1458 fse3t_b(:,:,:) = fse3t_n(:,:,:) 846 1459 neuler = 0 847 ELSE IF( id 2> 0 ) THEN1460 ELSE IF( id(2) > 0 ) THEN 848 1461 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_b not found in restart files' 849 1462 IF(lwp) write(numout,*) 'fse3t_b set equal to fse3t_n.' … … 867 1480 IF( ln_vvl_zstar ) THEN ! z_star case ! 868 1481 ! ! ----------- ! 869 IF( MIN( id 3, id4) > 0 ) THEN1482 IF( MIN( id(3), id(4) ) > 0 ) THEN 870 1483 CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' ) 871 1484 ENDIF … … 873 1486 ELSE ! z_tilde and layer cases ! 874 1487 ! ! ----------------------- ! 875 IF( MIN( id 3, id4) > 0 ) THEN ! all required arrays exist1488 IF( MIN( id(3), id(4) ) > 0 ) THEN ! all required arrays exist 876 1489 CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:) ) 877 1490 CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:) ) … … 879 1492 tilde_e3t_b(:,:,:) = 0.0_wp 880 1493 tilde_e3t_n(:,:,:) = 0.0_wp 881 ENDIF 882 ! ! ------------ ! 1494 ! 1495 neuler = 0 1496 ENDIF ! ------------ ! 883 1497 IF( ln_vvl_ztilde ) THEN ! z_tilde case ! 884 1498 ! ! ------------ ! 885 IF( id5 > 0 ) THEN ! required array exists 886 CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:) ) 887 ELSE ! array is missing 888 hdiv_lf(:,:,:) = 0.0_wp 1499 IF( MINVAL(id(5:7) ) > 0 ) THEN ! all required arrays exist 1500 CALL iom_get( numror, jpdom_autoglo, 'hdivn_lf', hdivn_lf(:,:,:) ) 1501 CALL iom_get( numror, jpdom_autoglo, 'un_lf' , un_lf(:,:,:) ) 1502 CALL iom_get( numror, jpdom_autoglo, 'vn_lf' , vn_lf(:,:,:) ) 1503 ELSE ! one at least array is missing 1504 hdivn_lf(:,:,:) = 0.0_wp 1505 un_lf(:,:,:) = 0.0_wp 1506 vn_lf(:,:,:) = 0.0_wp 1507 neuler = 0 889 1508 ENDIF 890 1509 ENDIF 1510 ! 891 1511 ENDIF 892 1512 ! … … 898 1518 tilde_e3t_b(:,:,:) = 0.0_wp 899 1519 tilde_e3t_n(:,:,:) = 0.0_wp 900 IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0.0_wp901 1520 END IF 1521 IF( ln_vvl_ztilde ) THEN 1522 hdivn_lf(:,:,:) = 0.0_wp 1523 un_lf(:,:,:) = 0.0_wp 1524 vn_lf(:,:,:) = 0.0_wp 1525 ENDIF 902 1526 ENDIF 903 1527 … … 919 1543 IF( ln_vvl_ztilde ) THEN ! z_tilde case ! 920 1544 ! ! ------------ ! 921 CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:) ) 1545 CALL iom_rstput( kt, nitrst, numrow, 'hdivn_lf', hdivn_lf(:,:,:) ) 1546 CALL iom_rstput( kt, nitrst, numrow, 'un_lf' , un_lf(:,:,:) ) 1547 CALL iom_rstput( kt, nitrst, numrow, 'vn_lf' , vn_lf(:,:,:) ) 922 1548 ENDIF 923 1549 … … 926 1552 927 1553 END SUBROUTINE dom_vvl_rst 928 929 1554 930 1555 SUBROUTINE dom_vvl_ctl … … 938 1563 INTEGER :: ios 939 1564 940 NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, & 941 & ln_vvl_zstar_at_eqtor , rn_ahe3 , rn_rst_e3t , & 942 & rn_lf_cutoff , rn_zdef_max , ln_vvl_dbg ! not yet implemented: ln_vvl_kepe 1565 NAMELIST/nam_vvl/ ln_vvl_zstar , ln_vvl_ztilde , & 1566 & ln_vvl_layer , ln_vvl_ztilde_as_zstar , & 1567 & ln_vvl_zstar_at_eqtor , ln_vvl_zstar_on_shelf , & 1568 & ln_vvl_adv_cn2 , ln_vvl_adv_fct , & 1569 & ln_vvl_lap , ln_vvl_blp , & 1570 & rn_ahe3_lap , rn_ahe3_blp , & 1571 & rn_rst_e3t , rn_lf_cutoff , & 1572 & ln_vvl_regrid , rn_zdef_max , & 1573 & ln_vvl_dbg ! not yet implemented: ln_vvl_kepe 943 1574 !!---------------------------------------------------------------------- 944 1575 … … 961 1592 WRITE(numout,*) ' layer ln_vvl_layer = ', ln_vvl_layer 962 1593 WRITE(numout,*) ' ztilde as zstar ln_vvl_ztilde_as_zstar = ', ln_vvl_ztilde_as_zstar 963 WRITE(numout,*) ' ztilde near the equator ln_vvl_zstar_at_eqtor = ', ln_vvl_zstar_at_eqtor964 1594 ! WRITE(numout,*) ' Namelist nam_vvl : chose kinetic-to-potential energy conservation' 965 1595 ! WRITE(numout,*) ' ln_vvl_kepe = ', ln_vvl_kepe 966 WRITE(numout,*) ' Namelist nam_vvl : thickness diffusion coefficient' 967 WRITE(numout,*) ' rn_ahe3 = ', rn_ahe3 1596 WRITE(numout,*) ' ztilde near the equator ln_vvl_zstar_at_eqtor = ', ln_vvl_zstar_at_eqtor 1597 WRITE(numout,*) ' ztilde on shelves ln_vvl_zstar_on_shelf = ', ln_vvl_zstar_on_shelf 1598 WRITE(numout,*) ' Namelist nam_vvl : thickness advection scheme' 1599 WRITE(numout,*) ' 2nd order ln_vvl_adv_cn2 = ', ln_vvl_adv_cn2 1600 WRITE(numout,*) ' 2nd order FCT ln_vvl_adv_fct = ', ln_vvl_adv_fct 1601 WRITE(numout,*) ' Namelist nam_vvl : thickness diffusion scheme' 1602 WRITE(numout,*) ' Laplacian ln_vvl_lap = ', ln_vvl_lap 1603 WRITE(numout,*) ' Bilaplacian ln_vvl_blp = ', ln_vvl_blp 1604 WRITE(numout,*) ' Laplacian coefficient rn_ahe3_lap = ', rn_ahe3_lap 1605 WRITE(numout,*) ' Bilaplacian coefficient rn_ahe3_blp = ', rn_ahe3_blp 1606 WRITE(numout,*) ' Namelist nam_vvl : layers regriding' 1607 WRITE(numout,*) ' ln_vvl_regrid = ', ln_vvl_regrid 968 1608 WRITE(numout,*) ' Namelist nam_vvl : maximum e3t deformation fractional change' 969 1609 WRITE(numout,*) ' rn_zdef_max = ', rn_zdef_max … … 992 1632 993 1633 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' ) 994 IF( .NOT. ln_vvl_zstar .AND. nn_isf .NE. 0) CALL ctl_stop( 'Only vvl_zstar has been tested with ice shelf cavity' ) 1634 1635 IF ( ln_vvl_ztilde.OR.ln_vvl_layer ) THEN 1636 ioptio = 0 ! Choose one advection scheme at most 1637 IF( ln_vvl_adv_cn2 ) ioptio = ioptio + 1 1638 IF( ln_vvl_adv_fct ) ioptio = ioptio + 1 1639 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE thickness advection scheme in namelist nam_vvl' ) 1640 ENDIF 995 1641 996 1642 IF(lwp) THEN ! Print the choice … … 1401 2047 END SUBROUTINE dom_vvl_orca_fix 1402 2048 2049 SUBROUTINE dom_vvl_zdf( kt, p2dt ) 2050 !!---------------------------------------------------------------------- 2051 !! *** ROUTINE dom_vvl_zdf *** 2052 !! 2053 !! ** Purpose : Do vertical thicknesses anomaly diffusion 2054 !! 2055 !! ** Method : 2056 !! 2057 !! ** Action : 2058 !!--------------------------------------------------------------------- 2059 USE oce , ONLY: zwd => ua , zws => va ! (ua,va) used as 3D workspace 2060 ! 2061 INTEGER , INTENT(in) :: kt ! ocean time-step index 2062 REAL(wp) , INTENT(in) :: p2dt ! time step 2063 ! 2064 INTEGER :: ji, jj, jk ! dummy loop indices 2065 REAL(wp) :: zr_tscale 2066 REAL(wp) :: za1, za2, za3, za4, zdiff 2067 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwt 2068 !!--------------------------------------------------------------------- 2069 ! 2070 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_zdf') 2071 ! 2072 zr_tscale = 0.5 2073 ! 2074 CALL wrk_alloc( jpi, jpj, jpk, zwi, zwt) 2075 ! 2076 ! 2077 2078 zwt(:,:,1:jpk) = 1.e-10 2079 2080 zwt(:,:,1) = 0.0_wp 2081 ! DO jk = 2, jpkm1 2082 ! DO jj = 1, jpj 2083 ! DO ji = 1, jpi 2084 ! zwt(ji,jj,jk) = zwt(ji,jj,jk-1) + (tilde_e3t_a(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1) 2085 ! END DO 2086 ! END DO 2087 ! END DO 2088 ! 2089 ! 2090 ! Set diffusivity (homogeneous to an inverse time scale) 2091 ! 2092 DO jk = 2, jpkm1 2093 DO jj = 2, jpjm1 2094 DO ji = fs_2, fs_jpim1 2095 ! Taper a little bit: 2096 za1 = tilde_e3t_n(ji,jj,jk-1)+e3t_0(ji,jj,jk-1) 2097 za2 = tilde_e3t_n(ji,jj,jk )+e3t_0(ji,jj,jk ) 2098 za4 = 0.5_wp * (e3t_0(ji,jj,jk-1) + e3t_0(ji,jj,jk )) 2099 za3 = 0.5_wp * (za1 + za2) 2100 zdiff = ABS(za3-za4)/za4 2101 ! IF (zdiff>=0.8) THEN 2102 ! zwt(ji,jj,jk) = zr_tscale * MIN(zdiff,1._wp) * za3 / p2dt * tmask(ji,jj,jk) 2103 zwt(ji,jj,jk) = dsm(ji,jj)/ht_0(ji,jj)*(1._wp-tanh((mbkt(ji,jj)+1-jk)*0.2))*tmask(ji,jj,jk) 2104 2105 ! ELSE 2106 ! zwt(ji,jj,jk) = 0.e0*tmask(ji,jj,jk) 2107 ! ENDIF 2108 END DO 2109 END DO 2110 END DO 2111 ! 2112 ! 2113 ! Diagonal, lower (i), upper (s) (including the bottom boundary condition since avt is masked) 2114 DO jk = 1, jpkm1 2115 DO jj = 2, jpjm1 2116 DO ji = fs_2, fs_jpim1 ! vector opt. 2117 zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk ) / fse3w(ji,jj,jk ) 2118 zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / fse3w(ji,jj,jk+1) 2119 zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk) 2120 END DO 2121 END DO 2122 END DO 2123 ! 2124 !! Matrix inversion from the first level 2125 !!---------------------------------------------------------------------- 2126 DO jj = 2, jpjm1 2127 DO ji = fs_2, fs_jpim1 2128 zwt(ji,jj,1) = zwd(ji,jj,1) 2129 END DO 2130 END DO 2131 DO jk = 2, jpkm1 2132 DO jj = 2, jpjm1 2133 DO ji = fs_2, fs_jpim1 2134 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1) 2135 END DO 2136 END DO 2137 END DO 2138 ! 2139 ! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1 2140 DO jk = 2, jpkm1 2141 DO jj = 2, jpjm1 2142 DO ji = fs_2, fs_jpim1 2143 tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * tilde_e3t_a(ji,jj,jk-1) 2144 END DO 2145 END DO 2146 END DO 2147 ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 2148 DO jj = 2, jpjm1 2149 DO ji = fs_2, fs_jpim1 2150 tilde_e3t_a(ji,jj,jpkm1) = tilde_e3t_a(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 2151 END DO 2152 END DO 2153 DO jk = jpk-2, 1, -1 2154 DO jj = 2, jpjm1 2155 DO ji = fs_2, fs_jpim1 2156 tilde_e3t_a(ji,jj,jk) = ( tilde_e3t_a(ji,jj,jk) - zws(ji,jj,jk) * tilde_e3t_a(ji,jj,jk+1) ) & 2157 & / zwt(ji,jj,jk) * tmask(ji,jj,jk) 2158 END DO 2159 END DO 2160 END DO 2161 ! 2162 CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwt) 2163 ! 2164 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_zdf') 2165 ! 2166 END SUBROUTINE dom_vvl_zdf 2167 2168 SUBROUTINE dom_vvl_zdf2( kt, p2dt ) 2169 !!---------------------------------------------------------------------- 2170 !! *** ROUTINE dom_vvl_zdf *** 2171 !! 2172 !! ** Purpose : Do vertical interface diffusion 2173 !! 2174 !! ** Method : 2175 !! 2176 !! ** Action : 2177 !!--------------------------------------------------------------------- 2178 USE oce , ONLY: zwd => ua , zws => va ! (ua,va) used as 3D workspace 2179 ! 2180 INTEGER , INTENT(in) :: kt ! ocean time-step index 2181 REAL(wp) , INTENT(in) :: p2dt ! time step 2182 ! 2183 INTEGER :: ji, jj, jk, kbot, kbotm1 ! dummy loop indices 2184 REAL(wp) :: zr_tscale 2185 REAL(wp) :: za1, za2, za3, za4, zdiff 2186 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwt, zw 2187 !!--------------------------------------------------------------------- 2188 ! 2189 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_zdf') 2190 ! 2191 zr_tscale = 0.5 2192 ! 2193 CALL wrk_alloc( jpi, jpj, jpk, zwi, zwt, zw) 2194 ! 2195 ! Compute internal interfaces depths: 2196 zw(:,:,1) = 0._wp 2197 DO jk = 2, jpk 2198 DO jj = 2, jpjm1 2199 DO ji = fs_2, fs_jpim1 ! vector opt. 2200 zw(ji,jj,jk) = zw(ji,jj,jk-1) + (tilde_e3t_a(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1) 2201 END DO 2202 END DO 2203 END DO 2204 ! 2205 ! Set diffusivities at interfaces 2206 zwt(:,:,:) = 0.00000001_wp * tmask(:,:,:) 2207 zwt(:,:,1) = 0._wp 2208 ! 2209 ! 2210 ! Diagonal, lower (i), upper (s) (including the bottom boundary condition since avt is masked) 2211 DO jk = 2, jpkm1 2212 DO jj = 2, jpjm1 2213 DO ji = fs_2, fs_jpim1 ! vector opt. 2214 zwi(ji,jj,jk) = - 0.5_wp * p2dt * (zwt(ji,jj,jk-1) + zwt(ji,jj,jk )) & 2215 & / fse3w(ji,jj,jk-1) / fse3t(ji,jj,jk-1) & 2216 & * ht(ji,jj) 2217 zws(ji,jj,jk) = - 0.5_wp * p2dt * (zwt(ji,jj,jk ) + zwt(ji,jj,jk+1)) & 2218 & / fse3w(ji,jj,jk ) / fse3t(ji,jj,jk ) & 2219 & * ht(ji,jj) 2220 2221 zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk) 2222 END DO 2223 END DO 2224 END DO 2225 ! Boundary conditions (Neumann) 2226 DO jj = 2, jpjm1 2227 DO ji = fs_2, fs_jpim1 2228 zwi(ji,jj,1) = 0._wp 2229 zws(ji,jj,1) = 0._wp 2230 zwd(ji,jj,1) = 1._wp 2231 zw (ji,jj,1) = 0._wp 2232 ! 2233 zwd(ji,jj,2) = zwd(ji,jj,2) + zwi(ji,jj,2) 2234 zwi(ji,jj,2) = 0._wp 2235 ! zwi(ji,jj,2) = 0._wp 2236 ! zws(ji,jj,2) = 0._wp 2237 ! zwd(ji,jj,2) = 1._wp 2238 END DO 2239 END DO 2240 DO jj = 2, jpjm1 2241 DO ji = fs_2, fs_jpim1 2242 kbot = mbkt(ji,jj) + 1 2243 kbotm1 = mbkt(ji,jj) 2244 zwi(ji,jj,kbot ) = 0._wp 2245 zws(ji,jj,kbot ) = 0._wp 2246 zwd(ji,jj,kbot ) = 1._wp 2247 ! 2248 zwd(ji,jj,kbotm1) = zwd(ji,jj,kbotm1) + zws(ji,jj,kbotm1) 2249 zws(ji,jj,kbotm1) = 0._wp 2250 ! zwi(ji,jj,kbotm1) = 0._wp 2251 ! zws(ji,jj,kbotm1) = 0._wp 2252 ! zwd(ji,jj,kbotm1) = 1._wp 2253 END DO 2254 END DO 2255 ! 2256 DO jk = 2, jpkm1 2257 DO jj = 2, jpjm1 2258 DO ji = fs_2, fs_jpim1 2259 zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 2260 END DO 2261 END DO 2262 END DO 2263 ! 2264 DO jk = 2, jpk 2265 DO jj = 2, jpjm1 2266 DO ji = fs_2, fs_jpim1 ! vector opt. 2267 zwi(ji,jj,jk) = zw(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) *zwi(ji,jj,jk-1) 2268 END DO 2269 END DO 2270 END DO 2271 ! 2272 DO jk = jpk-1, 2, -1 2273 DO jj = 2, jpjm1 2274 DO ji = fs_2, fs_jpim1 ! vector opt. 2275 zw(ji,jj,jk) = ( zwi(ji,jj,jk) - zws(ji,jj,jk) * zw(ji,jj,jk+1) ) / zwd(ji,jj,jk) 2276 END DO 2277 END DO 2278 END DO 2279 ! 2280 ! Revert to thicknesses anomalies: 2281 DO jk = 1, jpkm1 2282 DO jj = 2, jpjm1 2283 DO ji = fs_2, fs_jpim1 ! vector opt. 2284 tilde_e3t_a(ji,jj,jk) = (zw(ji,jj,jk+1)-zw(ji,jj,jk)-e3t_0(ji,jj,jk))* tmask(ji,jj,jk) 2285 END DO 2286 END DO 2287 END DO 2288 ! 2289 CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwt, zw) 2290 ! 2291 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_zdf') 2292 ! 2293 END SUBROUTINE dom_vvl_zdf2 2294 2295 SUBROUTINE dom_vvl_regrid( kt ) 2296 !!---------------------------------------------------------------------- 2297 !! *** ROUTINE dom_vvl_regrid *** 2298 !! 2299 !! ** Purpose : Ensure "well-behaved" vertical grid 2300 !! 2301 !! ** Method : More or less adapted from references below. 2302 !! 2303 !! ** Action : Ensure that thickness are above a given value, spaced enough 2304 !! and revert to Eulerian coordinates near the bottom. 2305 !! 2306 !! References : Bleck, R. and S. Benjamin, 1993: Regional Weather Prediction 2307 !! with a Model Combining Terrain-following and Isentropic 2308 !! coordinates. Part I: Model Description. Monthly Weather Rev., 2309 !! 121, 1770-1785. 2310 !! Toy, M., 2011: Incorporating Condensational Heating into a 2311 !! Nonhydrostatic Atmospheric Model Based on a Hybrid Isentropic- 2312 !! Sigma Vertical Coordinate. Monthly Weather Rev., 139, 2940-2954. 2313 !!---------------------------------------------------------------------- 2314 !! * Arguments 2315 INTEGER, INTENT( in ) :: kt ! time step 2316 2317 !! * Local declarations 2318 INTEGER :: ji, jj, jk ! dummy loop indices 2319 LOGICAL :: ll_chk_bot2top, ll_chk_top2bot, ll_lapdiff_cond 2320 LOGICAL :: ll_e3tsurf_const, ll_zdiff_cond, ll_blpdiff_cond 2321 INTEGER :: jkbot 2322 REAL(wp) :: zh_min, zh_0, zh2, zdiff, zh_max, ztmph, ztmpd, dzmin_vvl 2323 REAL(wp) :: zufim1, zufi, zvfjm1, zvfj 2324 REAL(wp) :: zh_new, zh_old, zh_bef, ztmp, ztmp1, z2dt, zh_up, zh_dwn 2325 REAL(wp) :: zeu2, zev2, zfrch_stp, zfrch_rel, zfrac_bot, zscal_bot 2326 REAL(wp) :: zhdiff, zhdiff2, zvdiff, zhlim, zhlim2, zvlim 2327 REAL(wp), POINTER, DIMENSION(:,:) :: zdw 2328 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwdw, zwdw_b 2329 !!---------------------------------------------------------------------- 2330 2331 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_regrid') 2332 ! 2333 CALL wrk_alloc( jpi, jpj, jpk, zwdw) 2334 ! 2335 ! Some user defined parameters below: 2336 ll_chk_bot2top = .TRUE. 2337 ll_chk_top2bot = .TRUE. 2338 ll_e3tsurf_const = .FALSE. 2339 dzmin_vvl = 1._wp ! Absolute minimum depth (in meters) 2340 zfrch_stp = 0.5_wp ! Maximum fractionnal thickness change in one time step (<= 1.) 2341 zfrch_rel = 0.0_wp ! Maximum relative thickness change in the vertical (<= 1.) 2342 zfrac_bot = 0.01_wp ! Fraction of bottom level allowed to change 2343 zscal_bot = 2._wp ! Depth lengthscale 2344 ll_zdiff_cond = .FALSE. ! Conditionnal vertical diffusion of interfaces 2345 zvdiff = 0.1_wp ! m/s 2346 zvlim = 0.4_wp ! max d2h/dh 2347 ll_lapdiff_cond = .FALSE. ! Conditionnal Laplacian diffusion of interfaces 2348 zhdiff = 0.1_wp ! ad. 2349 zhlim = 0.01_wp ! max lap(z)/e1 2350 ll_blpdiff_cond = .FALSE. ! Conditionnal Bilaplacian diffusion of interfaces 2351 zhdiff2 = 0.8_wp ! ad. 2352 zhlim2 = 5.e-11_wp ! max bilap(z)/e1**3 2353 ! --------------------------------------------------------------------------------------- 2354 ! 2355 ! Set arrays determining maximum vertical displacement at the bottom: 2356 !-------------------------------------------------------------------- 2357 IF ( kt==nit000 ) THEN 2358 DO jj = 2, jpjm1 2359 DO ji = 2, jpim1 2360 jk = MIN(mbkt(ji,jj), mbkt(ji+1,jj), mbkt(ji-1,jj), mbkt(ji,jj+1), mbkt(ji,jj-1)) 2361 i_int_bot(ji,jj) = jk 2362 END DO 2363 END DO 2364 dsm(:,:) = REAL( i_int_bot(:,:), wp ) ; CALL lbc_lnk(dsm(:,:),'T',1.) 2365 i_int_bot(:,:) = MAX( INT( dsm(:,:) ), 1 ) 2366 2367 CALL wrk_alloc( jpi, jpj, zdw ) 2368 DO jj = 2, jpjm1 2369 DO ji = 2, jpim1 2370 zdw(ji,jj) = MAX(ABS(ht_0(ji,jj)-ht_0(ji+1,jj))*umask(ji ,jj,1), & 2371 & ABS(ht_0(ji,jj)-ht_0(ji-1,jj))*umask(ji-1,jj,1), & 2372 & ABS(ht_0(ji,jj)-ht_0(ji,jj+1))*vmask(ji,jj ,1), & 2373 & ABS(ht_0(ji,jj)-ht_0(ji,jj-1))*vmask(ji,jj-1,1) ) 2374 zdw(ji,jj) = MAX(zscal_bot * zdw(ji,jj), rsmall ) 2375 END DO 2376 END DO 2377 CALL lbc_lnk( zdw(:,:), 'T', 1. ) 2378 2379 DO jj = 2, jpjm1 2380 DO ji = 2, jpim1 2381 dsm(ji,jj) = 1._wp/16._wp * ( zdw(ji-1,jj-1) + zdw(ji+1,jj-1) & 2382 & + zdw(ji-1,jj+1) + zdw(ji+1,jj+1) & 2383 & + 2._wp*( zdw(ji ,jj-1) + zdw(ji-1,jj ) & 2384 & + zdw(ji+1,jj ) + zdw(ji ,jj+1) ) & 2385 & + 4._wp* zdw(ji ,jj ) ) 2386 END DO 2387 END DO 2388 2389 CALL lbc_lnk( dsm(:,:), 'T', 1. ) 2390 CALL wrk_dealloc( jpi, jpj, zdw) 2391 2392 IF (ln_zps) THEN 2393 DO jj = 1, jpj 2394 DO ji = 1, jpi 2395 jk = i_int_bot(ji,jj) 2396 hsm(ji,jj) = zfrac_bot * e3w_1d(jk) 2397 END DO 2398 END DO 2399 ELSE 2400 DO jj = 1, jpj 2401 DO ji = 1, jpi 2402 jk = i_int_bot(ji,jj) 2403 hsm(ji,jj) = zfrac_bot * e3w_0(ji,jj,jk) 2404 END DO 2405 END DO 2406 ENDIF 2407 END IF 2408 2409 ! Provisionnal interface depths: 2410 !------------------------------- 2411 zwdw(:,:,1) = 0.e0 2412 DO jj = 1, jpj 2413 DO ji = 1, jpi 2414 DO jk = 2, jpk 2415 zwdw(ji,jj,jk) = zwdw(ji,jj,jk-1) + & 2416 & (tilde_e3t_a(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1) 2417 END DO 2418 END DO 2419 END DO 2420 ! 2421 ! Conditionnal horizontal Laplacian diffusion: 2422 !--------------------------------------------- 2423 IF ( ll_lapdiff_cond ) THEN 2424 CALL wrk_alloc( jpi, jpj, jpk, zwdw_b) 2425 ! 2426 zwdw_b(:,:,1) = 0._wp 2427 DO jj = 1, jpj 2428 DO ji = 1, jpi 2429 DO jk=2,jpk 2430 zwdw_b(ji,jj,jk) = zwdw_b(ji,jj,jk-1) + & 2431 & (tilde_e3t_b(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1) 2432 END DO 2433 END DO 2434 END DO 2435 ! 2436 DO jk = 2, jpkm1 2437 DO jj = 1, jpjm1 2438 DO ji = 1, fs_jpim1 ! vector opt. 2439 ua(ji,jj,jk) = umask(ji,jj,jk) * re2u_e1u(ji,jj) & 2440 & * ( zwdw_b(ji,jj,jk) - zwdw_b(ji+1,jj ,jk) ) 2441 va(ji,jj,jk) = vmask(ji,jj,jk) * re1v_e2v(ji,jj) & 2442 & * ( zwdw_b(ji,jj,jk) - zwdw_b(ji ,jj+1,jk) ) 2443 END DO 2444 END DO 2445 END DO 2446 2447 DO jk = 2, jpkm1 2448 DO jj = 2, jpjm1 2449 DO ji = fs_2, fs_jpim1 ! vector opt. 2450 ztmp1 = ( (ua(ji-1,jj ,jk) - ua(ji,jj,jk)) & 2451 & + (va(ji ,jj-1,jk) - va(ji,jj,jk)) ) * r1_e12t(ji,jj) 2452 zh2 = MAX(abs(ztmp1)-zhlim*SQRT(r1_e12t(ji,jj)), 0._wp) 2453 ztmp = SIGN(zh2, ztmp1) 2454 zeu2 = zhdiff * e12t(ji,jj)*e12t(ji,jj)/(e1t(ji,jj)*e1t(ji,jj) + e2t(ji,jj)*e2t(ji,jj)) 2455 zwdw(ji,jj,jk) = zwdw(ji,jj,jk) + zeu2 * ztmp * tmask(ji,jj,jk) 2456 END DO 2457 END DO 2458 END DO 2459 ! 2460 CALL wrk_dealloc( jpi, jpj, jpk, zwdw_b) 2461 ! 2462 ENDIF 2463 2464 ! Conditionnal horizontal Bilaplacian diffusion: 2465 !----------------------------------------------- 2466 IF ( ll_blpdiff_cond ) THEN 2467 CALL wrk_alloc( jpi, jpj, jpk, zwdw_b) 2468 ! 2469 zwdw_b(:,:,1) = 0._wp 2470 DO jj = 1, jpj 2471 DO ji = 1, jpi 2472 DO jk=2,jpk 2473 zwdw_b(ji,jj,jk) = zwdw_b(ji,jj,jk-1) + & 2474 & (tilde_e3t_b(ji,jj,jk-1)+e3t_0(ji,jj,jk-1)) * tmask(ji,jj,jk-1) 2475 END DO 2476 END DO 2477 END DO 2478 ! 2479 DO jk = 2, jpkm1 2480 DO jj = 1, jpjm1 2481 DO ji = 1, fs_jpim1 ! vector opt. 2482 ua(ji,jj,jk) = umask(ji,jj,jk) * re2u_e1u(ji,jj) & 2483 & * ( zwdw_b(ji,jj,jk) - zwdw_b(ji+1,jj ,jk) ) 2484 va(ji,jj,jk) = vmask(ji,jj,jk) * re1v_e2v(ji,jj) & 2485 & * ( zwdw_b(ji,jj,jk) - zwdw_b(ji ,jj+1,jk) ) 2486 END DO 2487 END DO 2488 END DO 2489 2490 DO jk = 2, jpkm1 2491 DO jj = 2, jpjm1 2492 DO ji = fs_2, fs_jpim1 ! vector opt. 2493 zwdw_b(ji,jj,jk) = -( (ua(ji-1,jj ,jk) - ua(ji,jj,jk)) & 2494 & + (va(ji ,jj-1,jk) - va(ji,jj,jk)) ) * r1_e12t(ji,jj) 2495 END DO 2496 END DO 2497 END DO 2498 ! 2499 CALL lbc_lnk( zwdw_b(:,:,:), 'T', 1. ) 2500 ! 2501 DO jk = 2, jpkm1 2502 DO jj = 1, jpjm1 2503 DO ji = 1, fs_jpim1 ! vector opt. 2504 ua(ji,jj,jk) = umask(ji,jj,jk) * re2u_e1u(ji,jj) & 2505 & * ( zwdw_b(ji,jj,jk) - zwdw_b(ji+1,jj ,jk) ) 2506 va(ji,jj,jk) = vmask(ji,jj,jk) * re1v_e2v(ji,jj) & 2507 & * ( zwdw_b(ji,jj,jk) - zwdw_b(ji ,jj+1,jk) ) 2508 END DO 2509 END DO 2510 END DO 2511 ! 2512 DO jk = 2, jpkm1 2513 DO jj = 2, jpjm1 2514 DO ji = fs_2, fs_jpim1 ! vector opt. 2515 ztmp1 = ( (ua(ji-1,jj ,jk) - ua(ji,jj,jk)) & 2516 & + (va(ji ,jj-1,jk) - va(ji,jj,jk)) ) * r1_e12t(ji,jj) 2517 zh2 = MAX(abs(ztmp1)-zhlim2, 0._wp) 2518 ztmp = SIGN(zh2, ztmp1) 2519 zeu2 = zhdiff2 * e12t(ji,jj)*e12t(ji,jj) / 16._wp 2520 zwdw(ji,jj,jk) = zwdw(ji,jj,jk) + zeu2 * ztmp * tmask(ji,jj,jk) 2521 END DO 2522 END DO 2523 END DO 2524 ! 2525 CALL wrk_dealloc( jpi, jpj, jpk, zwdw_b) 2526 ! 2527 ENDIF 2528 2529 ! Conditionnal vertical diffusion: 2530 !--------------------------------- 2531 IF ( ll_zdiff_cond ) THEN 2532 DO jk = 2, jpkm1 2533 DO jj = 2, jpjm1 2534 DO ji = fs_2, fs_jpim1 ! vector opt. 2535 ztmp = -( (tilde_e3t_b(ji,jj,jk-1)+e3t_0(ji,jj,jk-1))*tmask(ji,jj,jk-1) & 2536 -(tilde_e3t_b(ji,jj,jk )+e3t_0(ji,jj,jk ))*tmask(ji,jj,jk ) ) 2537 ztmp1 = 0.5_wp * ( tilde_e3t_b(ji,jj,jk-1) + e3t_0(ji,jj,jk-1) & 2538 & +tilde_e3t_b(ji,jj,jk ) + e3t_0(ji,jj,jk ) ) 2539 zh2 = MAX(abs(ztmp)-zvlim*ztmp1, 0._wp) 2540 ztmp = SIGN(zh2, ztmp) 2541 IF ((jk==mbkt(ji,jj)).AND.(ln_zps)) ztmp=0.e0 2542 zwdw(ji,jj,jk) = zwdw(ji,jj,jk) + zvdiff * ztmp * tmask(ji,jj,jk) 2543 END DO 2544 END DO 2545 END DO 2546 ENDIF 2547 ! 2548 ! Check grid from the bottom to the surface 2549 !------------------------------------------ 2550 IF ( ll_chk_bot2top ) THEN 2551 DO jj = 2, jpjm1 2552 DO ji = 2, jpim1 2553 jkbot = mbkt(ji,jj) 2554 DO jk = jkbot,2,-1 2555 ! 2556 zh_0 = e3t_0(ji,jj,jk) 2557 zh_bef = MIN(tilde_e3t_b(ji,jj,jk) + zh_0, tilde_e3t_b(ji,jj,jk-1) + e3t_0(ji,jj,jk-1)) 2558 zh_old = tilde_e3t_a(ji,jj,jk ) + zh_0 2559 ! zh_dwn = tilde_e3t_a(ji,jj,jk+1) + e3t_0(ji,jj,jk+1) 2560 zh_min = MIN(zh_0/3._wp, dzmin_vvl) 2561 ! 2562 ! Set maximum and minimum vertical excursions 2563 ztmph = hsm(ji,jj) 2564 ztmpd = dsm(ji,jj) 2565 zh2 = ztmph * exp(-(gdepw_0(ji,jj,jk)-gdepw_0(ji,jj,i_int_bot(ji,jj)+1))/ztmpd) 2566 zdiff = cush_max(gdepw_0(ji,jj,jk)-zwdw(ji,jj,jk), zh2 ) 2567 zwdw(ji,jj,jk) = MAX(zwdw(ji,jj,jk), gdepw_0(ji,jj,jk) - zdiff) 2568 zdiff = cush_max(zwdw(ji,jj,jk)-gdepw_0(ji,jj,jk), zh2 ) 2569 zwdw(ji,jj,jk) = MIN(zwdw(ji,jj,jk), gdepw_0(ji,jj,jk) + zdiff) 2570 ! 2571 ! New layer thickness: 2572 zh_new = zwdw(ji,jj,jk+1) - zwdw(ji,jj,jk) 2573 ! 2574 ! Ensure minimum layer thickness: 2575 ! zh_new = MIN(zh_old, zh_dwn * zfrch_rel / (2._wp-zfrch_rel) ) 2576 zh_new = cush(zh_new, zh_min) 2577 ! 2578 ! Final flux: 2579 zdiff = (zh_new - zh_old) * tmask(ji,jj,jk) 2580 ! 2581 ! Limit flux: 2582 ztmp = MIN( ABS(zdiff), zfrch_stp*zh_bef ) 2583 zdiff = SIGN(ztmp, zh_new - zh_old) 2584 zh_new = zdiff + zh_old 2585 ! 2586 tilde_e3t_a(ji,jj,jk ) = (zh_new - e3t_0(ji,jj,jk)) * tmask(ji,jj,jk) 2587 zwdw(ji,jj,jk) = zwdw(ji,jj,jk+1) - zh_new 2588 tilde_e3t_a(ji,jj,jk-1) = (-zdiff + tilde_e3t_a(ji,jj,jk-1) ) * tmask(ji,jj,jk-1) 2589 END DO 2590 END DO 2591 END DO 2592 END IF 2593 ! 2594 ! Check grid from the surface to the bottom 2595 !------------------------------------------ 2596 IF ( ll_chk_top2bot ) THEN 2597 DO jj = 2, jpjm1 2598 DO ji = 2, jpim1 2599 jkbot = mbkt(ji,jj) 2600 DO jk = 1, jkbot-1 2601 ! 2602 zh_0 = e3t_0(ji,jj,jk) 2603 zh_bef = MIN(tilde_e3t_b(ji,jj,jk) + zh_0, tilde_e3t_b(ji,jj,jk+1) + e3t_0(ji,jj,jk+1)) 2604 zh_old = tilde_e3t_a(ji,jj,jk ) + zh_0 2605 ! zh_up = tilde_e3t_a(ji,jj,jk-1) + e3t_0(ji,jj,jk-1) 2606 zh_min = MIN(zh_0/3._wp, dzmin_vvl) 2607 IF ((jk<=5).AND.ll_e3tsurf_const) zh_min = MAX(e3t_0(ji,jj,1)/3._wp, dzmin_vvl) 2608 ! 2609 ! Ensure minimum layer thickness: 2610 ! zh_new=MIN(zh_old, zh_up * zfrch_rel / (2._wp-zfrch_rel) ) 2611 zh_new = cush(zh_old, zh_min) 2612 ! 2613 ! Final flux: 2614 zdiff = (zh_new -zh_old) * tmask(ji,jj,jk) 2615 ! 2616 ! Limit flux: 2617 ztmp = MIN( ABS(zdiff), zfrch_stp*zh_bef ) 2618 zdiff = SIGN(ztmp, zdiff) 2619 zh_new = zdiff + zh_old 2620 ! 2621 tilde_e3t_a(ji,jj,jk ) = (zh_new - e3t_0(ji,jj,jk)) * tmask(ji,jj,jk) 2622 zwdw(ji,jj,jk+1) = zwdw(ji,jj,jk) + zh_new 2623 tilde_e3t_a(ji,jj,jk+1) = (-zdiff + tilde_e3t_a(ji,jj,jk+1) ) * tmask(ji,jj,jk+1) 2624 END DO 2625 ! 2626 END DO 2627 END DO 2628 ENDIF 2629 ! 2630 CALL wrk_dealloc( jpi, jpj, jpk, zwdw ) 2631 ! 2632 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_regrid') 2633 ! 2634 END SUBROUTINE dom_vvl_regrid 2635 2636 FUNCTION cush(hin, hmin) RESULT(hout) 2637 !!---------------------------------------------------------------------- 2638 !! *** FUNCTION cush *** 2639 !! 2640 !! ** Purpose : 2641 !! 2642 !! ** Method : 2643 !! 2644 !!---------------------------------------------------------------------- 2645 IMPLICIT NONE 2646 REAL(wp), INTENT(in) :: hin, hmin 2647 REAL(wp) :: hout, zx, zh_cri 2648 !!---------------------------------------------------------------------- 2649 zh_cri = 3._wp * hmin 2650 ! 2651 IF ( hin<=0._wp ) THEN 2652 hout = hmin 2653 ! 2654 ELSEIF ( (hin>0._wp).AND.(hin<=zh_cri) ) THEN 2655 zx = hin/zh_cri 2656 hout = hmin * (1._wp + zx + zx*zx) 2657 ! 2658 ELSEIF ( hin>zh_cri ) THEN 2659 hout = hin 2660 ! 2661 ENDIF 2662 ! 2663 END FUNCTION cush 2664 2665 FUNCTION cush_max(hin, hmax) RESULT(hout) 2666 !!---------------------------------------------------------------------- 2667 !! *** FUNCTION cush *** 2668 !! 2669 !! ** Purpose : 2670 !! 2671 !! ** Method : 2672 !! 2673 !!---------------------------------------------------------------------- 2674 IMPLICIT NONE 2675 REAL(wp), INTENT(in) :: hin, hmax 2676 REAL(wp) :: hout, hmin, zx, zh_cri 2677 !!---------------------------------------------------------------------- 2678 hmin = 0.1_wp * hmax 2679 zh_cri = 3._wp * hmin 2680 ! 2681 IF ( (hin>=(hmax-zh_cri)).AND.(hin<=(hmax-hmin))) THEN 2682 zx = (hmax-hin)/zh_cri 2683 hout = hmax - hmin * (1._wp + zx + zx*zx) 2684 ! 2685 ELSEIF ( hin>(hmax-zh_cri) ) THEN 2686 hout = hmax - hmin 2687 ! 2688 ELSE 2689 hout = hin 2690 ! 2691 ENDIF 2692 ! 2693 END FUNCTION cush_max 2694 2695 SUBROUTINE dom_vvl_adv_fct( kt, pta, uin, vin ) 2696 !!---------------------------------------------------------------------- 2697 !! *** ROUTINE dom_vvl_adv_fct *** 2698 !! 2699 !! ** Purpose : Do thickness advection 2700 !! 2701 !! ** Method : FCT scheme to ensure positivity 2702 !! 2703 !! ** Action : - Update pta thickness tendency and diffusive fluxes 2704 !! - this is the total trend, hence it does include sea level motions 2705 !! - Upstream corrections to antidiffusive fluxes ensure 2706 !! that barotropic transport matches what is contained in input fluxes 2707 !!---------------------------------------------------------------------- 2708 ! 2709 INTEGER, INTENT( in ) :: kt ! ocean time-step index 2710 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pta ! thickness baroclinic trend 2711 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in) :: uin, vin ! input velocities 2712 ! 2713 INTEGER :: ji, jj, jk ! dummy loop indices 2714 INTEGER :: ikbu, ikbv, ibot 2715 REAL(wp) :: z2dtt, zbtr, ztra ! local scalar 2716 REAL(wp) :: zdi, zdj, zmin ! - - 2717 REAL(wp) :: zfp_ui, zfp_vj ! - - 2718 REAL(wp) :: zfm_ui, zfm_vj ! - - 2719 REAL(wp) :: zfp_hi, zfp_hj ! - - 2720 REAL(wp) :: zfm_hi, zfm_hj ! - - 2721 REAL(wp), POINTER, DIMENSION(:,:) :: zbu, zbv, zhu_b, zhv_b 2722 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwx, zwy, zwi 2723 !!---------------------------------------------------------------------- 2724 ! 2725 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_adv_fct') 2726 ! 2727 CALL wrk_alloc( jpi, jpj, zhu_b, zhv_b, zbu, zbv) 2728 CALL wrk_alloc( jpi, jpj, jpk, zwx, zwy, zwi) 2729 ! 2730 ! 1. Initializations 2731 ! ------------------ 2732 ! 2733 IF( neuler == 0 .AND. kt == nit000 ) THEN 2734 z2dtt = rdt 2735 ELSE 2736 z2dtt = 2.0_wp * rdt 2737 ENDIF 2738 ! 2739 zwi(:,:,:) = 0.e0 2740 zwx(:,:,:) = 0.e0 2741 zwy(:,:,:) = 0.e0 2742 ! 2743 ! 2744 ! 2. upstream advection with initial mass fluxes & intermediate update 2745 ! -------------------------------------------------------------------- 2746 DO jk = 1, jpkm1 2747 DO jj = 1, jpjm1 2748 DO ji = 1, fs_jpim1 ! vector opt. 2749 ! 2750 zfp_ui = uin(ji,jj,jk) + ABS( uin(ji,jj,jk) ) 2751 zfm_ui = uin(ji,jj,jk) - ABS( uin(ji,jj,jk) ) 2752 zfp_hi = fse3t_b(ji ,jj ,jk) 2753 zfm_hi = fse3t_b(ji+1,jj ,jk) 2754 zwx(ji,jj,jk) = 0.5 * e2u(ji,jj) * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk) 2755 ! 2756 zfp_vj = vin(ji,jj,jk) + ABS( vin(ji,jj,jk) ) 2757 zfm_vj = vin(ji,jj,jk) - ABS( vin(ji,jj,jk) ) 2758 zfp_hj = fse3t_b(ji ,jj ,jk) 2759 zfm_hj = fse3t_b(ji ,jj+1,jk) 2760 zwy(ji,jj,jk) = 0.5 * e1v(ji,jj) * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk) 2761 END DO 2762 END DO 2763 END DO 2764 2765 IF ( .NOT.ln_sco ) THEN 2766 ! Correct bottom upstream fluxes 2767 ! considering "shelf horizon depths" 2768 DO jj = 1, jpjm1 2769 DO ji = 1, fs_jpim1 ! vector opt. 2770 ! upstream scheme 2771 ikbu = mbku(ji,jj) 2772 ikbv = mbkv(ji,jj) 2773 zfp_ui = uin(ji,jj,ikbu) + ABS( uin(ji,jj,ikbu) ) 2774 zfm_ui = uin(ji,jj,ikbu) - ABS( uin(ji,jj,ikbu) ) 2775 zfp_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji ,jj ,ikbu), 0._wp) 2776 zfm_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji+1,jj ,ikbu), 0._wp) 2777 zwx(ji,jj,ikbu) = 0.5 * e2u(ji,jj) * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,ikbu) 2778 ! 2779 zfp_vj = vin(ji,jj,ikbv) + ABS( vin(ji,jj,ikbv) ) 2780 zfm_vj = vin(ji,jj,ikbv) - ABS( vin(ji,jj,ikbv) ) 2781 zfp_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji ,jj ,ikbv), 0._wp) 2782 zfm_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji ,jj+1,ikbv), 0._wp) 2783 zwy(ji,jj,ikbv) = 0.5 * e1v(ji,jj) * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,ikbv) 2784 END DO 2785 END DO 2786 ENDIF 2787 2788 ! total advective trend 2789 DO jk = 1, jpkm1 2790 DO jj = 2, jpjm1 2791 DO ji = fs_2, fs_jpim1 ! vector opt. 2792 zbtr = r1_e12t(ji,jj) 2793 ! total intermediate advective trends 2794 ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 2795 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) 2796 ! 2797 ! update and guess with monotonic sheme 2798 pta(ji,jj,jk) = pta(ji,jj,jk) + ztra 2799 zwi(ji,jj,jk) = (fse3t_b(ji,jj,jk) + z2dtt * ztra ) * tmask(ji,jj,jk) 2800 END DO 2801 END DO 2802 END DO 2803 ! ! Lateral boundary conditions on zwi (unchanged sign) 2804 CALL lbc_lnk( zwi, 'T', 1. ) 2805 2806 IF ( ln_vvl_dbg ) THEN 2807 zmin = MINVAL( zwi(:,:,:), mask = tmask(:,:,:) == 1.e0 ) 2808 IF( lk_mpp ) CALL mpp_min( zmin ) 2809 IF( zmin < 0._wp) THEN 2810 IF(lwp) CALL ctl_stop('vvl_adv: CFL issue here') 2811 ENDIF 2812 ENDIF 2813 2814 ! 3. antidiffusive flux : high order minus low order 2815 ! -------------------------------------------------- 2816 ! antidiffusive flux on i and j 2817 DO jk = 1, jpkm1 2818 DO jj = 1, jpjm1 2819 DO ji = 1, fs_jpim1 ! vector opt. 2820 zwx(ji,jj,jk) = (e2u(ji,jj) * uin(ji,jj,jk) * fse3u_n(ji,jj,jk) & 2821 & - zwx(ji,jj,jk)) * umask(ji,jj,jk) 2822 zwy(ji,jj,jk) = (e1v(ji,jj) * vin(ji,jj,jk) * fse3v_n(ji,jj,jk) & 2823 & - zwy(ji,jj,jk)) * vmask(ji,jj,jk) 2824 ! 2825 ! Update advective fluxes 2826 un_td(ji,jj,jk) = un_td(ji,jj,jk) - zwx(ji,jj,jk) 2827 vn_td(ji,jj,jk) = vn_td(ji,jj,jk) - zwy(ji,jj,jk) 2828 END DO 2829 END DO 2830 END DO 2831 2832 CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) ! Lateral boundary conditions 2833 2834 ! 4. monotonicity algorithm 2835 ! ------------------------- 2836 CALL nonosc_2d( fse3t_b(:,:,:), zwx, zwy, zwi, z2dtt ) 2837 2838 ! 5. final trend with corrected fluxes 2839 ! ------------------------------------ 2840 ! 2841 ! Update advective fluxes 2842 un_td(:,:,:) = (un_td(:,:,:) + zwx(:,:,:))*umask(:,:,:) 2843 vn_td(:,:,:) = (vn_td(:,:,:) + zwy(:,:,:))*vmask(:,:,:) 2844 ! 2845 DO jk = 1, jpkm1 2846 DO jj = 2, jpjm1 2847 DO ji = fs_2, fs_jpim1 ! vector opt. 2848 ! 2849 zbtr = r1_e12t(ji,jj) 2850 ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 2851 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) 2852 ! add them to the general tracer trends 2853 pta(ji,jj,jk) = pta(ji,jj,jk) + ztra 2854 END DO 2855 END DO 2856 END DO 2857 2858 ! 2859 ! 6. Correct barotropic flux 2860 ! -------------------------- 2861 ! Compute barotropic flux difference: 2862 zbu(:,:) = 0.e0 2863 zbv(:,:) = 0.e0 2864 DO jj = 1, jpj 2865 DO ji = 1, jpi ! vector opt. 2866 DO jk = 1, jpkm1 2867 zbu(ji,jj) = zbu(ji,jj) - un_td(ji,jj,jk) * umask(ji,jj,jk) 2868 zbv(ji,jj) = zbv(ji,jj) - vn_td(ji,jj,jk) * vmask(ji,jj,jk) 2869 END DO 2870 END DO 2871 ENDDO 2872 2873 ! Compute upstream depths: 2874 zhu_b(:,:) = 0.e0 2875 zhv_b(:,:) = 0.e0 2876 IF ( ln_sco ) THEN; ibot=0 ; ELSE ; ibot=1 ; ENDIF 2877 2878 DO jj = 1, jpjm1 2879 DO ji = 1, jpim1 ! vector opt. 2880 zdi = 0.5_wp + 0.5_wp * SIGN(1._wp, zbu(ji,jj)) 2881 ikbu = mbku(ji,jj) - ibot 2882 DO jk = 1, ikbu 2883 zfp_hi = fse3t_b(ji ,jj ,jk) 2884 zfm_hi = fse3t_b(ji+1,jj ,jk) 2885 zhu_b(ji,jj) = zhu_b(ji,jj) + ( zdi * zfp_hi & 2886 & + (1._wp-zdi) * zfm_hi & 2887 & ) * umask(ji,jj,jk) 2888 END DO 2889 ! 2890 zdj = 0.5_wp + 0.5_wp * SIGN(1._wp, zbv(ji,jj)) 2891 ikbv = mbkv(ji,jj) - ibot 2892 DO jk = 1, ikbv 2893 zfp_hj = fse3t_b(ji ,jj ,jk) 2894 zfm_hj = fse3t_b(ji ,jj+1,jk) 2895 zhv_b(ji,jj) = zhv_b(ji,jj) + ( zdj * zfp_hj & 2896 & + (1._wp-zdj) * zfm_hj & 2897 & ) * vmask(ji,jj,jk) 2898 END DO 2899 END DO 2900 END DO 2901 2902 IF ( .NOT.ln_sco ) THEN 2903 ! Correct bottom value 2904 ! considering "shelf horizon depth" 2905 DO jj = 1, jpjm1 2906 DO ji = 1, jpim1 ! vector opt. 2907 zdi = 0.5_wp + 0.5_wp * SIGN(1._wp, zbu(ji,jj)) 2908 zdj = 0.5_wp + 0.5_wp * SIGN(1._wp, zbv(ji,jj)) 2909 ikbu = mbku(ji,jj) 2910 ikbv = mbkv(ji,jj) 2911 zfp_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji ,jj ,ikbu), 0._wp) 2912 zfm_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji+1,jj ,ikbu), 0._wp) 2913 zfp_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji ,jj ,ikbv), 0._wp) 2914 zfm_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji ,jj+1,ikbv), 0._wp) 2915 zhu_b(ji,jj) = zhu_b(ji,jj) + ( zdi * zfp_hi & 2916 & + (1._wp-zdi) * zfm_hi & 2917 & ) * umask(ji,jj,ikbu) 2918 zhv_b(ji,jj) = zhv_b(ji,jj) + ( zdj * zfp_hj & 2919 & + (1._wp-zdj) * zfm_hj & 2920 & ) * vmask(ji,jj,ikbv) 2921 END DO 2922 END DO 2923 ENDIF 2924 2925 ! Corrective barotropic velocity (times hor. scale factor) 2926 zbu(:,:) = zbu(:,:)/ (zhu_b(:,:)*umask(:,:,1)+1._wp-umask(:,:,1)) 2927 zbv(:,:) = zbv(:,:)/ (zhv_b(:,:)*vmask(:,:,1)+1._wp-vmask(:,:,1)) 2928 2929 CALL lbc_lnk( zbu(:,:), 'U', -1. ) 2930 CALL lbc_lnk( zbv(:,:), 'V', -1. ) 2931 CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) ! Lateral boundary conditions 2932 2933 ! Set corrective fluxes in upstream direction: 2934 ! 2935 zwx(:,:,:) = 0.e0 2936 zwy(:,:,:) = 0.e0 2937 DO jj = 1, jpjm1 2938 DO ji = 1, fs_jpim1 ! vector opt. 2939 ! upstream scheme 2940 zfp_ui = zbu(ji,jj) + ABS( zbu(ji,jj) ) 2941 zfm_ui = zbu(ji,jj) - ABS( zbu(ji,jj) ) 2942 zfp_vj = zbv(ji,jj) + ABS( zbv(ji,jj) ) 2943 zfm_vj = zbv(ji,jj) - ABS( zbv(ji,jj) ) 2944 DO jk = 1, jpkm1 2945 zfp_hi = fse3t_b(ji ,jj ,jk) 2946 zfm_hi = fse3t_b(ji+1,jj ,jk) 2947 ! 2948 zwx(ji,jj,jk) = 0.5 * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) * umask(ji,jj,jk) 2949 2950 zfp_hj = fse3t_b(ji ,jj ,jk) 2951 zfm_hj = fse3t_b(ji ,jj+1,jk) 2952 ! 2953 zwy(ji,jj,jk) = 0.5 * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) * vmask(ji,jj,jk) 2954 END DO 2955 END DO 2956 END DO 2957 2958 IF ( .NOT.ln_sco ) THEN 2959 ! Bottom correction: 2960 DO jj = 1, jpjm1 2961 DO ji = 1, fs_jpim1 ! vector opt. 2962 ! upstream scheme 2963 ikbu = mbku(ji,jj) 2964 ikbv = mbkv(ji,jj) 2965 ! 2966 zfp_ui = zbu(ji,jj) + ABS( zbu(ji,jj) ) 2967 zfm_ui = zbu(ji,jj) - ABS( zbu(ji,jj) ) 2968 zfp_vj = zbv(ji,jj) + ABS( zbv(ji,jj) ) 2969 zfm_vj = zbv(ji,jj) - ABS( zbv(ji,jj) ) 2970 ! 2971 zfp_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji ,jj ,ikbu), 0._wp) 2972 zfm_hi = MAX(hu_b(ji,jj) - fsdepw_b(ji+1,jj ,ikbu), 0._wp) 2973 ! 2974 zwx(ji,jj,ikbu) = 0.5 * ( zfp_ui * zfp_hi + zfm_ui * zfm_hi ) 2975 ! 2976 zfp_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji ,jj ,ikbv), 0._wp) 2977 zfm_hj = MAX(hv_b(ji,jj) - fsdepw_b(ji ,jj+1,ikbv), 0._wp) 2978 ! 2979 zwy(ji,jj,ikbv) = 0.5 * ( zfp_vj * zfp_hj + zfm_vj * zfm_hj ) 2980 END DO 2981 END DO 2982 ENDIF 2983 2984 CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) ! Lateral boundary conditions 2985 2986 un_td(:,:,:) = un_td(:,:,:) + zwx(:,:,:) 2987 vn_td(:,:,:) = vn_td(:,:,:) + zwy(:,:,:) 2988 ! 2989 ! Update trend with corrective fluxes: 2990 DO jk = 1, jpkm1 2991 DO jj = 2, jpjm1 2992 DO ji = fs_2, fs_jpim1 ! vector opt. 2993 ! 2994 zbtr = r1_e12t(ji,jj) 2995 ! total advective trends 2996 ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) & 2997 & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) 2998 ! add them to the general tracer trends 2999 pta(ji,jj,jk) = pta(ji,jj,jk) + ztra 3000 END DO 3001 END DO 3002 END DO 3003 3004 CALL wrk_dealloc( jpi, jpj, zhu_b, zhv_b, zbu, zbv) 3005 CALL wrk_dealloc( jpi, jpj, jpk, zwx, zwy, zwi) 3006 ! 3007 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_adv_fct') 3008 ! 3009 END SUBROUTINE dom_vvl_adv_fct 3010 3011 SUBROUTINE nonosc_2d( pbef, paa, pbb, paft, p2dt ) 3012 !!--------------------------------------------------------------------- 3013 !! *** ROUTINE nonosc_2d *** 3014 !! 3015 !! ** Purpose : compute monotonic thickness fluxes from the upstream 3016 !! scheme and the before field by a nonoscillatory algorithm 3017 !! 3018 !! ** Method : ... ??? 3019 !! warning : pbef and paft must be masked, but the boundaries 3020 !! conditions on the fluxes are not necessary zalezak (1979) 3021 !! drange (1995) multi-dimensional forward-in-time and upstream- 3022 !! in-space based differencing for fluid 3023 !!---------------------------------------------------------------------- 3024 ! 3025 !!---------------------------------------------------------------------- 3026 REAL(wp) , INTENT(in ) :: p2dt ! vertical profile of tracer time-step 3027 REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in ) :: pbef, paft ! before & after field 3028 REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) :: paa, pbb ! monotonic fluxes in the 3 directions 3029 ! 3030 INTEGER :: ji, jj, jk ! dummy loop indices 3031 REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt ! local scalars 3032 REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo ! - - 3033 REAL(wp) :: zupip1, zupim1, zupjp1, zupjm1, zupb, zupa 3034 REAL(wp) :: zdoip1, zdoim1, zdojp1, zdojm1, zdob, zdoa 3035 REAL(wp), POINTER, DIMENSION(:,:,:) :: zbetup, zbetdo, zbup, zbdo 3036 !!---------------------------------------------------------------------- 3037 ! 3038 IF( nn_timing == 1 ) CALL timing_start('nonosc2') 3039 ! 3040 CALL wrk_alloc( jpi, jpj, jpk, zbetup, zbetdo, zbup, zbdo ) 3041 ! 3042 3043 zbig = 1.e+40_wp 3044 zrtrn = 1.e-15_wp 3045 zbetup(:,:,jpk) = 0._wp ; zbetdo(:,:,jpk) = 0._wp 3046 3047 3048 ! Search local extrema 3049 ! -------------------- 3050 ! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land 3051 zbup = MAX( pbef * tmask - zbig * ( 1.e0 - tmask ), & 3052 & paft * tmask - zbig * ( 1.e0 - tmask ) ) 3053 zbdo = MIN( pbef * tmask + zbig * ( 1.e0 - tmask ), & 3054 & paft * tmask + zbig * ( 1.e0 - tmask ) ) 3055 3056 DO jk = 1, jpkm1 3057 z2dtt = p2dt 3058 DO jj = 2, jpjm1 3059 DO ji = fs_2, fs_jpim1 ! vector opt. 3060 3061 ! search maximum in neighbourhood 3062 zup = MAX( zbup(ji ,jj ,jk ), & 3063 & zbup(ji-1,jj ,jk ), zbup(ji+1,jj ,jk ), & 3064 & zbup(ji ,jj-1,jk ), zbup(ji ,jj+1,jk )) 3065 3066 ! search minimum in neighbourhood 3067 zdo = MIN( zbdo(ji ,jj ,jk ), & 3068 & zbdo(ji-1,jj ,jk ), zbdo(ji+1,jj ,jk ), & 3069 & zbdo(ji ,jj-1,jk ), zbdo(ji ,jj+1,jk )) 3070 3071 ! positive part of the flux 3072 zpos = MAX( 0., paa(ji-1,jj ,jk ) ) - MIN( 0., paa(ji ,jj ,jk ) ) & 3073 & + MAX( 0., pbb(ji ,jj-1,jk ) ) - MIN( 0., pbb(ji ,jj ,jk ) ) 3074 3075 ! negative part of the flux 3076 zneg = MAX( 0., paa(ji ,jj ,jk ) ) - MIN( 0., paa(ji-1,jj ,jk ) ) & 3077 & + MAX( 0., pbb(ji ,jj ,jk ) ) - MIN( 0., pbb(ji ,jj-1,jk ) ) 3078 3079 ! up & down beta terms 3080 zbt = e1t(ji,jj) * e2t(ji,jj) / z2dtt 3081 zbetup(ji,jj,jk) = ( zup - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt 3082 zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo ) / ( zneg + zrtrn ) * zbt 3083 END DO 3084 END DO 3085 END DO 3086 3087 CALL lbc_lnk( zbetup, 'T', 1. ) ; CALL lbc_lnk( zbetdo, 'T', 1. ) ! lateral boundary cond. (unchanged sign) 3088 3089 ! 3. monotonic flux in the i & j direction (paa & pbb) 3090 ! ---------------------------------------- 3091 DO jk = 1, jpkm1 3092 DO jj = 2, jpjm1 3093 DO ji = fs_2, fs_jpim1 ! vector opt. 3094 zau = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) ) 3095 zbu = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) ) 3096 zcu = ( 0.5 + SIGN( 0.5 , paa(ji,jj,jk) ) ) 3097 paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1.e0 - zcu) * zbu ) 3098 3099 zav = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) ) 3100 zbv = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) ) 3101 zcv = ( 0.5 + SIGN( 0.5 , pbb(ji,jj,jk) ) ) 3102 pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1.e0 - zcv) * zbv ) 3103 END DO 3104 END DO 3105 END DO 3106 CALL lbc_lnk( paa, 'U', -1. ) ; CALL lbc_lnk( pbb, 'V', -1. ) ! lateral boundary condition (changed sign) 3107 ! 3108 CALL wrk_dealloc( jpi, jpj, jpk, zbetup, zbetdo, zbup, zbdo ) 3109 ! 3110 IF( nn_timing == 1 ) CALL timing_stop('nonosc2') 3111 ! 3112 END SUBROUTINE nonosc_2d 3113 1403 3114 !!====================================================================== 1404 3115 END MODULE domvvl 1405 1406 1407
Note: See TracChangeset
for help on using the changeset viewer.