- Timestamp:
- 2020-08-21T18:26:57+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/OCE/stpMLF.F90
r13237 r13427 32 32 !! 4.0 ! 2017-05 (G. Madec) introduction of the vertical physics manager (zdfphy) 33 33 !! 4.1 ! 2019-08 (A. Coward, D. Storkey) rewrite in preparation for new timestepping scheme 34 !!---------------------------------------------------------------------- 35 36 !!---------------------------------------------------------------------- 37 !! stp_MLF : OPA system time-stepping 34 !! 4.x ! 2020-08 (S. Techene, G. Madec) quasi eulerian coordinate time stepping 35 !!---------------------------------------------------------------------- 36 37 #if defined key_qco 38 !!---------------------------------------------------------------------- 39 !! 'key_qco' Quasi-Eulerian vertical coordonate 40 !!---------------------------------------------------------------------- 41 42 !!---------------------------------------------------------------------- 43 !! stp_MLF : NEMO modified Leap Frog time-stepping with qco 38 44 !!---------------------------------------------------------------------- 39 45 USE step_oce ! time stepping definition modules … … 43 49 USE traatfqco ! time filtering (tra_atf_qco routine) 44 50 USE dynatfqco ! time filtering (dyn_atf_qco routine) 45 USE dynspg_ts ! surface pressure gradient: split-explicit scheme (define un_adv)51 46 52 USE bdydyn ! ocean open boundary conditions (define bdy_dyn) 47 53 … … 49 55 PRIVATE 50 56 51 #if ! defined key_qco52 !!----------------------------------------------------------------------53 !! 'key_qco' EMPTY MODULE Quasi-Eulerian vertical coordonate54 !!----------------------------------------------------------------------55 #else56 57 PUBLIC stp_MLF ! called by nemogcm.F90 57 58 … … 195 196 zgdept(:,:,jk) = gdept(:,:,jk,Nnn) 196 197 END DO 197 CALL ssh_nxt ( kstp, Nbb, Nnn, ssh, Naa ) ! after ssh (includes call to div_hor) 198 IF( .NOT.ln_linssh ) CALL dom_qco_r3c ( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa) ) ! "after" ssh./h._0 ratio 199 CALL wzv ( kstp, Nbb, Nnn, Naa, ww ) ! Nnn cross-level velocity 200 IF( ln_zad_Aimp ) CALL wAimp ( kstp, Nnn ) ! Adaptive-implicit vertical advection partitioning 201 CALL eos ( ts(:,:,:,:,Nnn), rhd, rhop, zgdept ) ! now in situ density for hpg computation 198 CALL ssh_nxt ( kstp, Nbb, Nnn, ssh, Naa ) ! after ssh (includes call to div_hor) 199 IF( .NOT.ln_linssh ) THEN 200 CALL dom_qco_r3c( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa) ) ! "after" ssh/h_0 ratio at t,u,v pts 201 IF( ln_dynspg_exp ) CALL dom_qco_r3c( ssh(:,:,Nnn), r3t(:,:,Nnn), r3u(:,:,Nnn), r3v(:,:,Nnn), r3f(:,:) ) ! spg_exp : needed only for "now" ssh/h_0 ratio at f point 202 ENDIF 203 CALL wzv ( kstp, Nbb, Nnn, Naa, ww ) ! Nnn cross-level velocity 204 IF( ln_zad_Aimp ) CALL wAimp ( kstp, Nnn ) ! Adaptive-implicit vertical advection partitioning 205 CALL eos ( ts(:,:,:,:,Nnn), rhd, rhop, zgdept ) ! now in situ density for hpg computation 202 206 203 207 … … 218 222 CALL dyn_hpg( kstp, Nnn , uu, vv, Nrhs ) ! horizontal gradient of Hydrostatic pressure 219 223 CALL dyn_spg( kstp, Nbb, Nnn, Nrhs, uu, vv, ssh, uu_b, vv_b, Naa ) ! surface pressure gradient 220 221 224 ! With split-explicit free surface, since now transports have been updated and ssh(:,:,Nrhs) as well 225 222 226 IF( ln_dynspg_ts ) THEN ! vertical scale factors and vertical velocity need to be updated 223 227 CALL div_hor ( kstp, Nbb, Nnn ) ! Horizontal divergence (2nd call in time-split case) 224 IF(.NOT.ln_linssh) CALL dom_qco_r3c ( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) ) ! "after" ssh./h._0 ratio228 IF(.NOT.ln_linssh) CALL dom_qco_r3c ( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) ) ! update ssh/h_0 ratio at t,u,v,f pts 225 229 ENDIF 226 230 CALL dyn_zdf ( kstp, Nbb, Nnn, Nrhs, uu, vv, Naa ) ! vertical diffusion … … 305 309 !! place. 306 310 !! 307 CALL mlf_baro_corr ( Nnn, Naa, uu, vv ) ! barotrope ajustment308 CALL finalize_ sbc ( kstp, Nbb, Naa, uu, vv, ts ) ! boundary condifions309 CALL tra_atf_qco ( kstp, Nbb, Nnn, Naa , ts )! time filtering of "now" tracer arrays310 CALL dyn_atf_qco ( kstp, Nbb, Nnn, Naa, uu, vv ) ! time filtering of "now" velocities and scale factors311 r3t(:,:,Nnn) = r3t_f(:,:) 311 IF( ln_dynspg_ts ) CALL mlf_baro_corr ( Nnn, Naa, uu, vv ) ! barotrope adjustment 312 CALL finalize_lbc ( kstp, Nbb , Naa, uu, vv, ts ) ! boundary conditions 313 CALL tra_atf_qco ( kstp, Nbb, Nnn, Naa , ts ) ! time filtering of "now" tracer arrays 314 CALL dyn_atf_qco ( kstp, Nbb, Nnn, Naa, uu, vv ) ! time filtering of "now" velocities 315 r3t(:,:,Nnn) = r3t_f(:,:) ! update now ssh/h_0 with time filtered values 312 316 r3u(:,:,Nnn) = r3u_f(:,:) 313 317 r3v(:,:,Nnn) = r3v_f(:,:) … … 347 351 ! Control 348 352 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 349 CALL stp_ctl ( kstp, N bb, Nnn, indic)353 CALL stp_ctl ( kstp, Nnn ) 350 354 351 355 IF( kstp == nit000 ) THEN ! 1st time step only … … 364 368 IF( kstp == nitend .OR. indic < 0 ) THEN 365 369 CALL iom_context_finalize( cxios_context ) ! needed for XIOS+AGRIF 366 IF( lrxios) CALL iom_context_finalize( crxios_context )370 IF( lrxios ) CALL iom_context_finalize( crxios_context ) 367 371 IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) ! 368 372 ENDIF … … 380 384 381 385 382 SUBROUTINE mlf_baro_corr (Kmm, Kaa, puu, pvv)386 SUBROUTINE mlf_baro_corr( Kmm, Kaa, puu, pvv ) 383 387 !!---------------------------------------------------------------------- 384 388 !! *** ROUTINE mlf_baro_corr *** … … 389 393 !! estimate (ln_dynspg_ts=T) 390 394 !! 391 !! ** Action : puu(Kmm),pvv(Kmm),puu(Kaa),pvv(Kaa) now and after horizontal velocity 392 !!---------------------------------------------------------------------- 393 INTEGER , INTENT(in ) :: Kmm, Kaa ! before and after time level indices 394 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! velocities 395 ! 396 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zue, zve 395 !! ** Action : puu(Kmm),pvv(Kmm) updated now horizontal velocity (ln_bt_fw=F) 396 !! puu(Kaa),pvv(Kaa) after horizontal velocity 397 !!---------------------------------------------------------------------- 398 USE dynspg_ts, ONLY : un_adv, vn_adv ! updated Kmm barotropic transport 399 !! 400 INTEGER , INTENT(in ) :: Kmm, Kaa ! before and after time level indices 401 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! velocities 397 402 ! 398 403 INTEGER :: jk ! dummy loop indices 399 !!---------------------------------------------------------------------- 400 401 IF ( ln_dynspg_ts ) THEN 402 ALLOCATE( zue(jpi,jpj) , zve(jpi,jpj) ) 403 ! Ensure below that barotropic velocities match time splitting estimate 404 ! Compute actual transport and replace it with ts estimate at "after" time step 405 zue(:,:) = e3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1) 406 zve(:,:) = e3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1) 407 DO jk = 2, jpkm1 408 zue(:,:) = zue(:,:) + e3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk) 409 zve(:,:) = zve(:,:) + e3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk) 404 REAL(wp), DIMENSION(jpi,jpj) :: zue, zve 405 !!---------------------------------------------------------------------- 406 407 ! Ensure below that barotropic velocities match time splitting estimate 408 ! Compute actual transport and replace it with ts estimate at "after" time step 409 zue(:,:) = e3u(:,:,1,Kaa) * puu(:,:,1,Kaa) * umask(:,:,1) 410 zve(:,:) = e3v(:,:,1,Kaa) * pvv(:,:,1,Kaa) * vmask(:,:,1) 411 DO jk = 2, jpkm1 412 zue(:,:) = zue(:,:) + e3u(:,:,jk,Kaa) * puu(:,:,jk,Kaa) * umask(:,:,jk) 413 zve(:,:) = zve(:,:) + e3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk) 414 END DO 415 DO jk = 1, jpkm1 416 puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - zue(:,:) * r1_hu(:,:,Kaa) + uu_b(:,:,Kaa) ) * umask(:,:,jk) 417 pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - zve(:,:) * r1_hv(:,:,Kaa) + vv_b(:,:,Kaa) ) * vmask(:,:,jk) 418 END DO 419 ! 420 IF( .NOT.ln_bt_fw ) THEN 421 ! Remove advective velocity from "now velocities" 422 ! prior to asselin filtering 423 ! In the forward case, this is done below after asselin filtering 424 ! so that asselin contribution is removed at the same time 425 DO jk = 1, jpkm1 426 puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm) + uu_b(:,:,Kmm) )*umask(:,:,jk) 427 pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm) + vv_b(:,:,Kmm) )*vmask(:,:,jk) 410 428 END DO 411 DO jk = 1, jpkm1412 puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - zue(:,:) * r1_hu(:,:,Kaa) + uu_b(:,:,Kaa) ) * umask(:,:,jk)413 pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - zve(:,:) * r1_hv(:,:,Kaa) + vv_b(:,:,Kaa) ) * vmask(:,:,jk)414 END DO415 !416 IF( .NOT.ln_bt_fw ) THEN417 ! Remove advective velocity from "now velocities"418 ! prior to asselin filtering419 ! In the forward case, this is done below after asselin filtering420 ! so that asselin contribution is removed at the same time421 DO jk = 1, jpkm1422 puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm) + uu_b(:,:,Kmm) )*umask(:,:,jk)423 pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm) + vv_b(:,:,Kmm) )*vmask(:,:,jk)424 END DO425 ENDIF426 !427 DEALLOCATE( zue, zve )428 !429 429 ENDIF 430 430 ! … … 432 432 433 433 434 SUBROUTINE finalize_ sbc (kt, Kbb, Kaa, puu, pvv, pts)435 !!---------------------------------------------------------------------- 436 !! *** ROUTINE finalize_ sbc ***434 SUBROUTINE finalize_lbc( kt, Kbb, Kaa, puu, pvv, pts ) 435 !!---------------------------------------------------------------------- 436 !! *** ROUTINE finalize_lbc *** 437 437 !! 438 438 !! ** Purpose : Apply the boundary condition on the after velocity … … 445 445 !! ** Action : puu(Kaa),pvv(Kaa) after horizontal velocity and tracers 446 446 !!---------------------------------------------------------------------- 447 INTEGER , INTENT(in ) :: kt ! ocean time-step index 448 INTEGER , INTENT(in ) :: Kbb, Kaa ! before and after time level indices 449 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! velocities to be time filtered 450 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers 447 INTEGER , INTENT(in ) :: kt ! ocean time-step index 448 INTEGER , INTENT(in ) :: Kbb, Kaa ! before and after time level indices 449 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt) , INTENT(inout) :: puu, pvv ! velocities to be time filtered 450 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers 451 !!---------------------------------------------------------------------- 451 452 ! 452 453 ! Update after tracer and velocity on domain lateral boundaries 453 454 ! 454 # if defined key_agrif455 CALL Agrif_tra ! AGRIF zoom boundaries456 CALL Agrif_dyn( kt ) !* AGRIF zoom boundaries457 # endif455 # if defined key_agrif 456 CALL Agrif_tra !* AGRIF zoom boundaries 457 CALL Agrif_dyn( kt ) 458 # endif 458 459 ! ! local domain boundaries (T-point, unchanged sign) 459 CALL lbc_lnk_multi( 'finalize_ sbc', puu(:,:,:, Kaa), 'U', -1., pvv(:,:,: ,Kaa), 'V', -1. &460 & , pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1. ) !* local domain boundaries460 CALL lbc_lnk_multi( 'finalize_lbc', puu(:,:,:, Kaa), 'U', -1., pvv(:,:,: ,Kaa), 'V', -1. & 461 & , pts(:,:,:,jp_tem,Kaa), 'T', 1., pts(:,:,:,jp_sal,Kaa), 'T', 1. ) 461 462 ! 462 463 ! !* BDY open boundaries … … 467 468 ENDIF 468 469 ! 469 END SUBROUTINE finalize_sbc 470 #endif 471 ! 470 END SUBROUTINE finalize_lbc 471 472 #else 473 !!---------------------------------------------------------------------- 474 !! default option EMPTY MODULE qco not activated 475 !!---------------------------------------------------------------------- 476 #endif 477 472 478 !!====================================================================== 473 479 END MODULE stepMLF
Note: See TracChangeset
for help on using the changeset viewer.
