Changeset 14789 for NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/OCE/DYN/dynspg_ts.F90
- Timestamp:
- 2021-05-05T13:18:04+02:00 (3 years ago)
- Location:
- NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/DYN/dynspg_ts.F90 (modified) (26 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev _r12970_AGRIF_CMEMSext/AGRIF5 ^/vendors/AGRIF/dev@HEAD ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 ^/vendors/PPR@HEAD ext/PPR 8 9 9 10 # SETTE 10 ^/utils/CI/sette@1 3559sette11 ^/utils/CI/sette@14244 sette
-
- Property svn:externals
-
NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/OCE/DYN/dynspg_ts.F90
r13546 r14789 117 117 118 118 119 SUBROUTINE dyn_spg_ts( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa )119 SUBROUTINE dyn_spg_ts( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa, k_only_ADV ) 120 120 !!---------------------------------------------------------------------- 121 121 !! … … 146 146 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 147 147 REAL(wp), DIMENSION(jpi,jpj,jpt) , INTENT(inout) :: pssh, puu_b, pvv_b ! SSH and barotropic velocities at main time levels 148 INTEGER , OPTIONAL , INTENT( in ) :: k_only_ADV ! only Advection in the RHS 148 149 ! 149 150 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 162 163 REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v ! top/bottom stress at u- & v-points 163 164 REAL(wp), DIMENSION(jpi,jpj) :: zhU, zhV ! fluxes 164 REAL(wp), DIMENSION(jpi, jpj, jpk) :: ze3u, ze3v 165 !!st#if defined key_qco 166 !!st REAL(wp), DIMENSION(jpi, jpj, jpk) :: ze3u, ze3v 167 !!st#endif 165 168 ! 166 169 REAL(wp) :: zwdramp ! local scalar - only used if ln_wd_dl = .True. … … 229 232 ! != zu_frc = 1/H e3*d/dt(Ua) =! (Vertical mean of Ua, the 3D trends) 230 233 ! ! --------------------------- ! 231 DO jk = 1 , jpk 232 ze3u(:,:,jk) = e3u(:,:,jk,Kmm)233 ze3v(:,:,jk) = e3v(:,:,jk,Kmm)234 END DO 235 !236 z u_frc(:,:) = SUM( ze3u(:,:,:) * uu(:,:,:,Krhs) * umask(:,:,:) , DIM=3 ) * r1_hu(:,:,Kmm)237 zv_frc(:,:) = SUM( ze3v(:,:,:) * vv(:,:,:,Krhs) * vmask(:,:,:) , DIM=3 ) * r1_hv(:,:,Kmm) 234 #if defined key_qco 235 zu_frc(:,:) = SUM( e3u_0(:,:,: ) * puu(:,:,:,Krhs) * umask(:,:,:), DIM=3 ) * r1_hu_0(:,:) 236 zv_frc(:,:) = SUM( e3v_0(:,:,: ) * pvv(:,:,:,Krhs) * vmask(:,:,:), DIM=3 ) * r1_hv_0(:,:) 237 #else 238 zu_frc(:,:) = SUM( e3u(:,:,:,Kmm) * puu(:,:,:,Krhs) * umask(:,:,:), DIM=3 ) * r1_hu(:,:,Kmm) 239 zv_frc(:,:) = SUM( e3v(:,:,:,Kmm) * pvv(:,:,:,Krhs) * vmask(:,:,:), DIM=3 ) * r1_hv(:,:,Kmm) 240 #endif 238 241 ! 239 242 ! 240 243 ! != U(Krhs) => baroclinic trend =! (remove its vertical mean) 241 244 DO jk = 1, jpkm1 ! ----------------------------- ! 242 uu(:,:,jk,Krhs) = (uu(:,:,jk,Krhs) - zu_frc(:,:) ) * umask(:,:,jk)243 vv(:,:,jk,Krhs) = (vv(:,:,jk,Krhs) - zv_frc(:,:) ) * vmask(:,:,jk)245 puu(:,:,jk,Krhs) = ( puu(:,:,jk,Krhs) - zu_frc(:,:) ) * umask(:,:,jk) 246 pvv(:,:,jk,Krhs) = ( pvv(:,:,jk,Krhs) - zv_frc(:,:) ) * vmask(:,:,jk) 244 247 END DO 245 248 … … 247 250 !!gm Is it correct to do so ? I think so... 248 251 249 ! != remove 2D Coriolis and pressure gradient trends=!250 ! ! -------------------------- -----------------------!252 ! != remove 2D Coriolis trend =! 253 ! ! -------------------------- ! 251 254 ! 252 255 IF( kt == nit000 .OR. .NOT. ln_linssh ) CALL dyn_cor_2D_init( Kmm ) ! Set zwz, the barotropic Coriolis force coefficient 253 ! ! recompute zwz = f/depth at every time step for (.NOT.ln_linssh) as the water colomn height changes 254 ! 255 ! !* 2D Coriolis trends 256 zhU(:,:) = puu_b(:,:,Kmm) * hu(:,:,Kmm) * e2u(:,:) ! now fluxes 257 zhV(:,:) = pvv_b(:,:,Kmm) * hv(:,:,Kmm) * e1v(:,:) ! NB: FULL domain : put a value in last row and column 258 ! 259 CALL dyn_cor_2d( ht(:,:), hu(:,:,Kmm), hv(:,:,Kmm), puu_b(:,:,Kmm), pvv_b(:,:,Kmm), zhU, zhV, & ! <<== in 260 & zu_trd, zv_trd ) ! ==>> out 261 ! 262 IF( .NOT.ln_linssh ) THEN !* surface pressure gradient (variable volume only) 263 ! 264 IF( ln_wd_il ) THEN ! W/D : limiter applied to spgspg 265 CALL wad_spg( pssh(:,:,Kmm), zcpx, zcpy ) ! Calculating W/D gravity filters, zcpx and zcpy 266 DO_2D( 0, 0, 0, 0 ) ! SPG with the application of W/D gravity filters 267 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( pssh(ji+1,jj ,Kmm) - pssh(ji ,jj ,Kmm) ) & 268 & * r1_e1u(ji,jj) * zcpx(ji,jj) * wdrampu(ji,jj) !jth 269 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( pssh(ji ,jj+1,Kmm) - pssh(ji ,jj ,Kmm) ) & 270 & * r1_e2v(ji,jj) * zcpy(ji,jj) * wdrampv(ji,jj) !jth 271 END_2D 272 ELSE ! now suface pressure gradient 273 DO_2D( 0, 0, 0, 0 ) 274 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( pssh(ji+1,jj ,Kmm) - pssh(ji ,jj ,Kmm) ) * r1_e1u(ji,jj) 275 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( pssh(ji ,jj+1,Kmm) - pssh(ji ,jj ,Kmm) ) * r1_e2v(ji,jj) 276 END_2D 277 ENDIF 278 ! 279 ENDIF 280 ! 281 DO_2D( 0, 0, 0, 0 ) ! Remove coriolis term (and possibly spg) from barotropic trend 282 zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj) 283 zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj) 284 END_2D 256 ! ! recompute zwz = f/depth at every time step for (.NOT.ln_linssh) as the water colomn height changes 257 ! 258 IF( .NOT. PRESENT(k_only_ADV) ) THEN !* remove the 2D Coriolis trend 259 zhU(:,:) = puu_b(:,:,Kmm) * hu(:,:,Kmm) * e2u(:,:) ! now fluxes 260 zhV(:,:) = pvv_b(:,:,Kmm) * hv(:,:,Kmm) * e1v(:,:) ! NB: FULL domain : put a value in last row and column 261 ! 262 CALL dyn_cor_2d( ht(:,:), hu(:,:,Kmm), hv(:,:,Kmm), puu_b(:,:,Kmm), pvv_b(:,:,Kmm), zhU, zhV, & ! <<== in 263 & zu_trd, zv_trd ) ! ==>> out 264 ! 265 DO_2D( 0, 0, 0, 0 ) ! Remove coriolis term (and possibly spg) from barotropic trend 266 zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj) 267 zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj) 268 END_2D 269 ENDIF 285 270 ! 286 271 ! != Add bottom stress contribution from baroclinic velocities =! 287 272 ! ! ----------------------------------------------------------- ! 288 CALL dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, zu_frc, zv_frc, zCdU_u, zCdU_v ) ! also provide the barotropic drag coefficients 273 IF( PRESENT(k_only_ADV) ) THEN !* only Advection in the RHS : provide the barotropic bottom drag coefficients 274 DO_2D( 0, 0, 0, 0 ) 275 zCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) 276 zCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) 277 END_2D 278 ELSE !* remove baroclinic drag AND provide the barotropic drag coefficients 279 CALL dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b, pvv_b, zu_frc, zv_frc, zCdU_u, zCdU_v ) 280 ENDIF 281 ! 289 282 ! != Add atmospheric pressure forcing =! 290 283 ! ! ---------------------------------- ! … … 306 299 ENDIF 307 300 ! 308 ! != Add atmospheric pressureforcing =!309 ! ! ------------------ ----------------!310 IF( ln_bt_fw ) THEN ! Add wind forcing301 ! != Add wind forcing =! 302 ! ! ------------------ ! 303 IF( ln_bt_fw ) THEN 311 304 DO_2D( 0, 0, 0, 0 ) 312 305 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rho0 * utau(ji,jj) * r1_hu(ji,jj,Kmm) … … 330 323 ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 fluxes (NOW) 331 324 zztmp = r1_rho0 * r1_2 332 zssh_frc(:,:) = zztmp * ( emp(:,:) + emp_b(:,:)&333 & - rnf(:,:) - rnf_b(:,:)&334 & + fwfisf_cav(:,:) + fwfisf_cav_b(:,:)&335 & + fwfisf_par(:,:) + fwfisf_par_b(:,:))325 zssh_frc(:,:) = zztmp * ( emp(:,:) + emp_b(:,:) & 326 & - rnf(:,:) - rnf_b(:,:) & 327 & + fwfisf_cav(:,:) + fwfisf_cav_b(:,:) & 328 & + fwfisf_par(:,:) + fwfisf_par_b(:,:) ) 336 329 ENDIF 337 330 ! != Add Stokes drift divergence =! (if exist) … … 386 379 ! 387 380 IF( ln_linssh ) THEN ! mid-step ocean depth is fixed (hup2_e=hu_n=hu_0) 388 zhup2_e(:,:) = hu (:,:,Kmm)389 zhvp2_e(:,:) = hv (:,:,Kmm)390 zhtp2_e(:,:) = ht (:,:)391 ENDIF 392 ! 393 IF (ln_bt_fw) THEN! FORWARD integration: start from NOW fields394 sshn_e(:,:) = pssh (:,:,Kmm)381 zhup2_e(:,:) = hu_0(:,:) 382 zhvp2_e(:,:) = hv_0(:,:) 383 zhtp2_e(:,:) = ht_0(:,:) 384 ENDIF 385 ! 386 IF( ln_bt_fw ) THEN ! FORWARD integration: start from NOW fields 387 sshn_e(:,:) = pssh (:,:,Kmm) 395 388 un_e (:,:) = puu_b(:,:,Kmm) 396 389 vn_e (:,:) = pvv_b(:,:,Kmm) … … 401 394 hvr_e (:,:) = r1_hv(:,:,Kmm) 402 395 ELSE ! CENTRED integration: start from BEFORE fields 403 sshn_e(:,:) = pssh (:,:,Kbb)396 sshn_e(:,:) = pssh (:,:,Kbb) 404 397 un_e (:,:) = puu_b(:,:,Kbb) 405 398 vn_e (:,:) = pvv_b(:,:,Kbb) … … 412 405 ! 413 406 ! Initialize sums: 414 puu_b 415 pvv_b 407 puu_b (:,:,Kaa) = 0._wp ! After barotropic velocities (or transport if flux form) 408 pvv_b (:,:,Kaa) = 0._wp 416 409 pssh (:,:,Kaa) = 0._wp ! Sum for after averaged sea level 417 un_adv(:,:) = 0._wp ! Sum for now transport issued from ts loop418 vn_adv(:,:) = 0._wp410 un_adv(:,:) = 0._wp ! Sum for now transport issued from ts loop 411 vn_adv(:,:) = 0._wp 419 412 ! 420 413 IF( ln_wd_dl ) THEN … … 475 468 ! 476 469 ! ! ocean u- and v-depth at mid-step (separate DO-loops remove the need of a lbc_lnk) 477 DO_2D( 1, 1, 1, 0 ) ! not jpi-column 470 #if defined key_qcoTest_FluxForm 471 ! ! 'key_qcoTest_FluxForm' : simple ssh average 472 DO_2D( 1, 0, 1, 1 ) ! not jpi-column 473 zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * ( zsshp2_e(ji,jj) + zsshp2_e(ji+1,jj ) ) * ssumask(ji,jj) 474 END_2D 475 DO_2D( 1, 1, 1, 0 ) 476 zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * ( zsshp2_e(ji,jj) + zsshp2_e(ji ,jj+1) ) * ssvmask(ji,jj) 477 END_2D 478 #else 479 ! ! no 'key_qcoTest_FluxForm' : surface weighted ssh average 480 DO_2D( 1, 0, 1, 1 ) ! not jpi-column 478 481 zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj) & 479 482 & * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 480 483 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj) 481 484 END_2D 482 DO_2D( 1, 0, 1, 1) ! not jpj-row485 DO_2D( 1, 1, 1, 0 ) ! not jpj-row 483 486 zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * r1_e1e2v(ji,jj) & 484 487 & * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 485 488 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj) 486 489 END_2D 490 #endif 487 491 ! 488 492 ENDIF … … 520 524 END_2D 521 525 ! 522 CALL lbc_lnk _multi( 'dynspg_ts', ssha_e, 'T', 1._wp, zhU, 'U', -1._wp, zhV, 'V', -1._wp )526 CALL lbc_lnk( 'dynspg_ts', ssha_e, 'T', 1._wp, zhU, 'U', -1._wp, zhV, 'V', -1._wp ) 523 527 ! 524 528 ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T) … … 540 544 ! 541 545 ! Sea Surface Height at u-,v-points (vvl case only) 542 IF( .NOT.ln_linssh ) THEN 546 IF( .NOT.ln_linssh ) THEN 547 #if defined key_qcoTest_FluxForm 548 ! ! 'key_qcoTest_FluxForm' : simple ssh average 549 DO_2D( 1, 0, 1, 1 ) 550 zsshu_a(ji,jj) = r1_2 * ( ssha_e(ji,jj) + ssha_e(ji+1,jj ) ) * ssumask(ji,jj) 551 END_2D 552 DO_2D( 1, 1, 1, 0 ) 553 zsshv_a(ji,jj) = r1_2 * ( ssha_e(ji,jj) + ssha_e(ji ,jj+1) ) * ssvmask(ji,jj) 554 END_2D 555 #else 543 556 DO_2D( 0, 0, 0, 0 ) 544 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 545 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 546 & + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) ) 547 zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 548 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 549 & + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) ) 550 END_2D 551 ENDIF 557 zsshu_a(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 558 & + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) ) * ssumask(ji,jj) 559 zsshv_a(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 560 & + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) ) * ssvmask(ji,jj) 561 END_2D 562 #endif 563 ENDIF 552 564 ! 553 565 ! Half-step back interpolation of SSH for surface pressure computation at step jit+1/2 … … 624 636 ! ! hu_e, hv_e hold depth at jn, zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2 625 637 ! ! backward interpolated depth used in spg terms at jn+1/2 638 #if defined key_qcoTest_FluxForm 639 ! ! 'key_qcoTest_FluxForm' : simple ssh average 640 zhu_bck = hu_0(ji,jj) + r1_2 * ( zsshp2_e(ji,jj) + zsshp2_e(ji+1,jj ) ) * ssumask(ji,jj) 641 zhv_bck = hv_0(ji,jj) + r1_2 * ( zsshp2_e(ji,jj) + zsshp2_e(ji ,jj+1) ) * ssvmask(ji,jj) 642 #else 626 643 zhu_bck = hu_0(ji,jj) + r1_2*r1_e1e2u(ji,jj) * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 627 644 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj) 628 645 zhv_bck = hv_0(ji,jj) + r1_2*r1_e1e2v(ji,jj) * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 629 646 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj) 647 #endif 630 648 ! ! inverse depth at jn+1 631 649 z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) ) … … 646 664 IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs 647 665 DO_2D( 0, 0, 0, 0 ) 648 ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rDt_e * zCdU_u(ji,jj) * hur_e(ji,jj))649 va_e(ji,jj) = va_e(ji,jj) /(1.0 - rDt_e * zCdU_v(ji,jj) * hvr_e(ji,jj))666 ua_e(ji,jj) = ua_e(ji,jj) / ( 1._wp - rDt_e * zCdU_u(ji,jj) * hur_e(ji,jj) ) 667 va_e(ji,jj) = va_e(ji,jj) / ( 1._wp - rDt_e * zCdU_v(ji,jj) * hvr_e(ji,jj) ) 650 668 END_2D 651 669 ENDIF 652 670 653 671 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 654 hu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1)655 hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / ( hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1))656 hv_e (2:jpim1,2:jpjm1) = hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1)657 hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / ( hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1))672 hu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1) 673 hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / ( hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1) ) 674 hv_e (2:jpim1,2:jpjm1) = hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1) 675 hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / ( hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1) ) 658 676 ENDIF 659 677 ! 660 678 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 661 CALL lbc_lnk _multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp &662 & 663 & 679 CALL lbc_lnk( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp & 680 & , hu_e , 'U', 1._wp, hv_e , 'V', 1._wp & 681 & , hur_e, 'U', 1._wp, hvr_e, 'V', 1._wp ) 664 682 ELSE 665 CALL lbc_lnk _multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp )683 CALL lbc_lnk( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp ) 666 684 ENDIF 667 685 ! ! open boundaries … … 743 761 ELSE 744 762 ! At this stage, pssh(:,:,:,Krhs) has been corrected: compute new depths at velocity points 763 #if defined key_qcoTest_FluxForm 764 ! ! 'key_qcoTest_FluxForm' : simple ssh average 745 765 DO_2D( 1, 0, 1, 0 ) 746 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 747 & * ( e1e2t(ji ,jj) * pssh(ji ,jj,Kaa) & 748 & + e1e2t(ji+1,jj) * pssh(ji+1,jj,Kaa) ) 749 zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 750 & * ( e1e2t(ji,jj ) * pssh(ji,jj ,Kaa) & 751 & + e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) ) 752 END_2D 753 CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 766 zsshu_a(ji,jj) = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji+1,jj ,Kaa) ) * ssumask(ji,jj) 767 zsshv_a(ji,jj) = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji ,jj+1,Kaa) ) * ssvmask(ji,jj) 768 END_2D 769 #else 770 DO_2D( 1, 0, 1, 0 ) 771 zsshu_a(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji ,jj) * pssh(ji ,jj,Kaa) & 772 & + e1e2t(ji+1,jj) * pssh(ji+1,jj,Kaa) ) * ssumask(ji,jj) 773 zsshv_a(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji,jj ) * pssh(ji,jj ,Kaa) & 774 & + e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) ) * ssvmask(ji,jj) 775 END_2D 776 #endif 777 CALL lbc_lnk( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 754 778 ! 755 779 DO jk=1,jpkm1 756 puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm) * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt_b 757 pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm) * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt_b 780 puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm) & 781 & * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt_b 782 pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm) & 783 & * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt_b 758 784 END DO 759 785 ! Save barotropic velocities not transport: … … 899 925 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 900 926 ! ! --------------- 901 IF( ln_rstart .AND. ln_bt_fw .AND. (.NOT.l_1st_euler)) THEN !* Read the restart file902 CALL iom_get( numror, jpdom_auto, 'ub2_b' , ub2_b (:,:), cd_type = 'U', psgn = -1._wp , ldxios = lrxios)903 CALL iom_get( numror, jpdom_auto, 'vb2_b' , vb2_b (:,:), cd_type = 'V', psgn = -1._wp , ldxios = lrxios)904 CALL iom_get( numror, jpdom_auto, 'un_bf' , un_bf (:,:), cd_type = 'U', psgn = -1._wp , ldxios = lrxios)905 CALL iom_get( numror, jpdom_auto, 'vn_bf' , vn_bf (:,:), cd_type = 'V', psgn = -1._wp , ldxios = lrxios)927 IF( ln_rstart .AND. ln_bt_fw .AND. .NOT.l_1st_euler ) THEN !* Read the restart file 928 CALL iom_get( numror, jpdom_auto, 'ub2_b' , ub2_b (:,:), cd_type = 'U', psgn = -1._wp ) 929 CALL iom_get( numror, jpdom_auto, 'vb2_b' , vb2_b (:,:), cd_type = 'V', psgn = -1._wp ) 930 CALL iom_get( numror, jpdom_auto, 'un_bf' , un_bf (:,:), cd_type = 'U', psgn = -1._wp ) 931 CALL iom_get( numror, jpdom_auto, 'vn_bf' , vn_bf (:,:), cd_type = 'V', psgn = -1._wp ) 906 932 IF( .NOT.ln_bt_av ) THEN 907 CALL iom_get( numror, jpdom_auto, 'sshbb_e' , sshbb_e(:,:), cd_type = 'T', psgn = 1._wp , ldxios = lrxios)908 CALL iom_get( numror, jpdom_auto, 'ubb_e' , ubb_e(:,:), cd_type = 'U', psgn = -1._wp , ldxios = lrxios)909 CALL iom_get( numror, jpdom_auto, 'vbb_e' , vbb_e(:,:), cd_type = 'V', psgn = -1._wp , ldxios = lrxios)910 CALL iom_get( numror, jpdom_auto, 'sshb_e' , sshb_e(:,:), cd_type = 'T', psgn = 1._wp , ldxios = lrxios)911 CALL iom_get( numror, jpdom_auto, 'ub_e' , ub_e(:,:), cd_type = 'U', psgn = -1._wp , ldxios = lrxios)912 CALL iom_get( numror, jpdom_auto, 'vb_e' , vb_e(:,:), cd_type = 'V', psgn = -1._wp , ldxios = lrxios)933 CALL iom_get( numror, jpdom_auto, 'sshbb_e' , sshbb_e(:,:), cd_type = 'T', psgn = 1._wp ) 934 CALL iom_get( numror, jpdom_auto, 'ubb_e' , ubb_e(:,:), cd_type = 'U', psgn = -1._wp ) 935 CALL iom_get( numror, jpdom_auto, 'vbb_e' , vbb_e(:,:), cd_type = 'V', psgn = -1._wp ) 936 CALL iom_get( numror, jpdom_auto, 'sshb_e' , sshb_e(:,:), cd_type = 'T', psgn = 1._wp ) 937 CALL iom_get( numror, jpdom_auto, 'ub_e' , ub_e(:,:), cd_type = 'U', psgn = -1._wp ) 938 CALL iom_get( numror, jpdom_auto, 'vb_e' , vb_e(:,:), cd_type = 'V', psgn = -1._wp ) 913 939 ENDIF 914 940 #if defined key_agrif 915 941 ! Read time integrated fluxes 916 942 IF ( .NOT.Agrif_Root() ) THEN 917 CALL iom_get( numror, jpdom_auto, 'ub2_i_b' , ub2_i_b(:,:), cd_type = 'U', psgn = -1._wp , ldxios = lrxios)918 CALL iom_get( numror, jpdom_auto, 'vb2_i_b' , vb2_i_b(:,:), cd_type = 'V', psgn = -1._wp , ldxios = lrxios)943 CALL iom_get( numror, jpdom_auto, 'ub2_i_b' , ub2_i_b(:,:), cd_type = 'U', psgn = -1._wp ) 944 CALL iom_get( numror, jpdom_auto, 'vb2_i_b' , vb2_i_b(:,:), cd_type = 'V', psgn = -1._wp ) 919 945 ELSE 920 946 ub2_i_b(:,:) = 0._wp ; vb2_i_b(:,:) = 0._wp ! used in the 1st update of agrif … … 935 961 ! ! ------------------- 936 962 IF(lwp) WRITE(numout,*) '---- ts_rst ----' 937 IF( lwxios ) CALL iom_swap( cwxios_context ) 938 CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:), ldxios = lwxios ) 939 CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:), ldxios = lwxios ) 940 CALL iom_rstput( kt, nitrst, numrow, 'un_bf' , un_bf (:,:), ldxios = lwxios ) 941 CALL iom_rstput( kt, nitrst, numrow, 'vn_bf' , vn_bf (:,:), ldxios = lwxios ) 963 CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:) ) 964 CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:) ) 965 CALL iom_rstput( kt, nitrst, numrow, 'un_bf' , un_bf (:,:) ) 966 CALL iom_rstput( kt, nitrst, numrow, 'vn_bf' , vn_bf (:,:) ) 942 967 ! 943 968 IF (.NOT.ln_bt_av) THEN 944 CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:) , ldxios = lwxios)945 CALL iom_rstput( kt, nitrst, numrow, 'ubb_e' , ubb_e(:,:) , ldxios = lwxios)946 CALL iom_rstput( kt, nitrst, numrow, 'vbb_e' , vbb_e(:,:) , ldxios = lwxios)947 CALL iom_rstput( kt, nitrst, numrow, 'sshb_e' , sshb_e(:,:) , ldxios = lwxios)948 CALL iom_rstput( kt, nitrst, numrow, 'ub_e' , ub_e(:,:) , ldxios = lwxios)949 CALL iom_rstput( kt, nitrst, numrow, 'vb_e' , vb_e(:,:) , ldxios = lwxios)969 CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:) ) 970 CALL iom_rstput( kt, nitrst, numrow, 'ubb_e' , ubb_e(:,:) ) 971 CALL iom_rstput( kt, nitrst, numrow, 'vbb_e' , vbb_e(:,:) ) 972 CALL iom_rstput( kt, nitrst, numrow, 'sshb_e' , sshb_e(:,:) ) 973 CALL iom_rstput( kt, nitrst, numrow, 'ub_e' , ub_e(:,:) ) 974 CALL iom_rstput( kt, nitrst, numrow, 'vb_e' , vb_e(:,:) ) 950 975 ENDIF 951 976 #if defined key_agrif 952 977 ! Save time integrated fluxes 953 978 IF ( .NOT.Agrif_Root() ) THEN 954 CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b' , ub2_i_b(:,:) , ldxios = lwxios)955 CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b' , vb2_i_b(:,:) , ldxios = lwxios)979 CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b' , ub2_i_b(:,:) ) 980 CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b' , vb2_i_b(:,:) ) 956 981 ENDIF 957 982 #endif 958 IF( lwxios ) CALL iom_swap( cxios_context )959 983 ENDIF 960 984 ! … … 1049 1073 CALL ts_rst( nit000, 'READ' ) 1050 1074 ! 1051 IF( lwxios ) THEN1052 ! define variables in restart file when writing with XIOS1053 CALL iom_set_rstw_var_active('ub2_b')1054 CALL iom_set_rstw_var_active('vb2_b')1055 CALL iom_set_rstw_var_active('un_bf')1056 CALL iom_set_rstw_var_active('vn_bf')1057 !1058 IF (.NOT.ln_bt_av) THEN1059 CALL iom_set_rstw_var_active('sshbb_e')1060 CALL iom_set_rstw_var_active('ubb_e')1061 CALL iom_set_rstw_var_active('vbb_e')1062 CALL iom_set_rstw_var_active('sshb_e')1063 CALL iom_set_rstw_var_active('ub_e')1064 CALL iom_set_rstw_var_active('vb_e')1065 ENDIF1066 #if defined key_agrif1067 ! Save time integrated fluxes1068 IF ( .NOT.Agrif_Root() ) THEN1069 CALL iom_set_rstw_var_active('ub2_i_b')1070 CALL iom_set_rstw_var_active('vb2_i_b')1071 ENDIF1072 #endif1073 ENDIF1074 !1075 1075 END SUBROUTINE dyn_spg_ts_init 1076 1076 … … 1086 1086 !! although they should be updated in the variable volume case. Not a big approximation. 1087 1087 !! To remove this approximation, copy lines below inside barotropic loop 1088 !! and update depths at T- F points (ht and zhf resp.) at each barotropic time step1088 !! and update depths at T- points (ht) at each barotropic time step 1089 1089 !! 1090 1090 !! Compute zwz = f / ( height of the water colomn ) … … 1093 1093 INTEGER :: ji ,jj, jk ! dummy loop indices 1094 1094 REAL(wp) :: z1_ht 1095 REAL(wp), DIMENSION(jpi,jpj) :: zhf1096 1095 !!---------------------------------------------------------------------- 1097 1096 ! 1098 1097 SELECT CASE( nvor_scheme ) 1099 CASE( np_EEN ) != EEN scheme using e3f energy & enstrophy scheme1100 SELECT CASE( nn_e en_e3f )!* ff_f/e3 at F-point1098 CASE( np_EEN, np_ENE, np_ENS , np_MIX ) != schemes using the same e3f definition 1099 SELECT CASE( nn_e3f_typ ) !* ff_f/e3 at F-point 1101 1100 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 1102 DO_2D( 1, 0, 1, 0 )1103 zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) +&1104 & ht(ji ,jj ) + ht(ji+1,jj )) * 0.25_wp1101 DO_2D( 0, 0, 0, 0 ) 1102 zwz(ji,jj) = ( ht(ji,jj+1) + ht(ji+1,jj+1) & 1103 & + ht(ji,jj ) + ht(ji+1,jj ) ) * 0.25_wp 1105 1104 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 1106 1105 END_2D 1107 1106 CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask) 1108 DO_2D( 1, 0, 1, 0 )1109 zwz(ji,jj) = ( ht (ji ,jj+1) + ht(ji+1,jj+1) &1110 & + ht (ji ,jj ) + ht(ji+1,jj ) ) &1111 & / ( MAX( 1._wp, ssmask(ji,jj+1) + ssmask(ji+1,jj+1) &1112 & + ssmask(ji ,jj ) + ssmask(ji+1,jj )) )1107 DO_2D( 0, 0, 0, 0 ) 1108 zwz(ji,jj) = ( ht(ji,jj+1) + ht(ji+1,jj+1) & 1109 & + ht(ji,jj ) + ht(ji+1,jj ) ) & 1110 & / ( MAX(ssmask(ji,jj+1) + ssmask(ji+1,jj+1) & 1111 & + ssmask(ji,jj ) + ssmask(ji+1,jj ) , 1._wp ) ) 1113 1112 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 1114 1113 END_2D 1115 1114 END SELECT 1116 1115 CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp ) 1117 ! 1118 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 1116 END SELECT 1117 ! 1118 SELECT CASE( nvor_scheme ) 1119 CASE( np_EEN ) 1120 ! 1121 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 1119 1122 DO_2D( 0, 1, 0, 1 ) 1120 1123 ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) … … 1124 1127 END_2D 1125 1128 ! 1126 CASE( np_EET ) != EEN scheme using e3t energy conserving scheme1127 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ;ftsw(1,:) = 0._wp1129 CASE( np_EET ) != EEN scheme using e3t energy conserving scheme 1130 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 1128 1131 DO_2D( 0, 1, 0, 1 ) 1129 1132 z1_ht = ssmask(ji,jj) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) ) … … 1134 1137 END_2D 1135 1138 ! 1136 CASE( np_ENE, np_ENS , np_MIX ) != all other schemes (ENE, ENS, MIX) except ENT !1137 !1138 zwz(:,:) = 0._wp1139 zhf(:,:) = 0._wp1140 1141 !!gm assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed1142 !!gm A priori a better value should be something like :1143 !!gm zhf(i,j) = masked sum of ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1)1144 !!gm divided by the sum of the corresponding mask1145 !!gm1146 !!1147 IF( .NOT.ln_sco ) THEN1148 1149 !!gm agree the JC comment : this should be done in a much clear way1150 1151 ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case1152 ! Set it to zero for the time being1153 ! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level1154 ! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth1155 ! ENDIF1156 ! zhf(:,:) = gdepw_0(:,:,jk+1)1157 !1158 ELSE1159 !1160 !zhf(:,:) = hbatf(:,:)1161 DO_2D( 1, 0, 1, 0 )1162 zhf(ji,jj) = ( ht_0 (ji,jj ) + ht_0 (ji+1,jj ) &1163 & + ht_0 (ji,jj+1) + ht_0 (ji+1,jj+1) ) &1164 & / MAX( ssmask(ji,jj ) + ssmask(ji+1,jj ) &1165 & + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp )1166 END_2D1167 ENDIF1168 !1169 DO jj = 1, jpjm11170 zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))1171 END DO1172 !1173 DO jk = 1, jpkm11174 DO jj = 1, jpjm11175 zhf(:,jj) = zhf(:,jj) + e3f(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)1176 END DO1177 END DO1178 CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp )1179 ! JC: TBC. hf should be greater than 01180 DO_2D( 1, 1, 1, 1 )1181 IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj)1182 END_2D1183 zwz(:,:) = ff_f(:,:) * zwz(:,:)1184 1139 END SELECT 1185 1140 1186 1141 END SUBROUTINE dyn_cor_2d_init 1187 1188 1142 1189 1143 … … 1308 1262 !!---------------------------------------------------------------------- 1309 1263 ! 1310 DO_2D( 1, 1, 1, 0) ! not jpi-column1264 DO_2D( 1, 0, 1, 1 ) ! not jpi-column 1311 1265 IF ( phU(ji,jj) > 0._wp ) THEN ; pUmsk(ji,jj) = pTmsk(ji ,jj) 1312 1266 ELSE ; pUmsk(ji,jj) = pTmsk(ji+1,jj) … … 1316 1270 END_2D 1317 1271 ! 1318 DO_2D( 1, 0, 1, 1) ! not jpj-row1272 DO_2D( 1, 1, 1, 0 ) ! not jpj-row 1319 1273 IF ( phV(ji,jj) > 0._wp ) THEN ; pVmsk(ji,jj) = pTmsk(ji,jj ) 1320 1274 ELSE ; pVmsk(ji,jj) = pTmsk(ji,jj+1) … … 1379 1333 END SUBROUTINE wad_spg 1380 1334 1381 1382 1335 1383 1336 SUBROUTINE dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, pu_RHSi, pv_RHSi, pCdU_u, pCdU_v )
Note: See TracChangeset
for help on using the changeset viewer.
