- Timestamp:
- 2020-09-29T12:41:06+02:00 (4 years ago)
- Location:
- NEMO/branches/2020/r12377_ticket2386
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/DYN/dynspg_ts.F90 (modified) (51 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/r12377_ticket2386
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEADext/AGRIF5 ^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 8 9 9 # SETTE 10 ^/utils/CI/sette@ HEADsette10 ^/utils/CI/sette@13507 sette
-
- Property svn:externals
-
NEMO/branches/2020/r12377_ticket2386/src/OCE/DYN/dynspg_ts.F90
r12511 r13540 87 87 !! * Substitutions 88 88 # include "do_loop_substitute.h90" 89 # include "domzgr_substitute.h90" 89 90 !!---------------------------------------------------------------------- 90 91 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 161 162 REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v ! top/bottom stress at u- & v-points 162 163 REAL(wp), DIMENSION(jpi,jpj) :: zhU, zhV ! fluxes 164 REAL(wp), DIMENSION(jpi, jpj, jpk) :: ze3u, ze3v 163 165 ! 164 166 REAL(wp) :: zwdramp ! local scalar - only used if ln_wd_dl = .True. … … 227 229 ! != zu_frc = 1/H e3*d/dt(Ua) =! (Vertical mean of Ua, the 3D trends) 228 230 ! ! --------------------------- ! 229 zu_frc(:,:) = SUM( e3u(:,:,:,Kmm) * uu(:,:,:,Krhs) * umask(:,:,:) , DIM=3 ) * r1_hu(:,:,Kmm) 230 zv_frc(:,:) = SUM( e3v(:,:,:,Kmm) * vv(:,:,:,Krhs) * vmask(:,:,:) , DIM=3 ) * r1_hv(:,:,Kmm) 231 DO jk = 1 , jpk 232 ze3u(:,:,jk) = e3u(:,:,jk,Kmm) 233 ze3v(:,:,jk) = e3v(:,:,jk,Kmm) 234 END DO 235 ! 236 zu_frc(:,:) = SUM( ze3u(:,:,:) * uu(:,:,:,Krhs) * umask(:,:,:) , DIM=3 ) * r1_hu(:,:,Kmm) 237 zv_frc(:,:) = SUM( ze3v(:,:,:) * vv(:,:,:,Krhs) * vmask(:,:,:) , DIM=3 ) * r1_hv(:,:,Kmm) 231 238 ! 232 239 ! … … 250 257 zhV(:,:) = pvv_b(:,:,Kmm) * hv(:,:,Kmm) * e1v(:,:) ! NB: FULL domain : put a value in last row and column 251 258 ! 252 CALL dyn_cor_2d( h u(:,:,Kmm), hv(:,:,Kmm), puu_b(:,:,Kmm), pvv_b(:,:,Kmm), zhU, zhV, & ! <<== in259 CALL dyn_cor_2d( ht(:,:), hu(:,:,Kmm), hv(:,:,Kmm), puu_b(:,:,Kmm), pvv_b(:,:,Kmm), zhU, zhV, & ! <<== in 253 260 & zu_trd, zv_trd ) ! ==>> out 254 261 ! … … 257 264 IF( ln_wd_il ) THEN ! W/D : limiter applied to spgspg 258 265 CALL wad_spg( pssh(:,:,Kmm), zcpx, zcpy ) ! Calculating W/D gravity filters, zcpx and zcpy 259 DO_2D _00_00266 DO_2D( 0, 0, 0, 0 ) ! SPG with the application of W/D gravity filters 260 267 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( pssh(ji+1,jj ,Kmm) - pssh(ji ,jj ,Kmm) ) & 261 268 & * r1_e1u(ji,jj) * zcpx(ji,jj) * wdrampu(ji,jj) !jth … … 264 271 END_2D 265 272 ELSE ! now suface pressure gradient 266 DO_2D _00_00273 DO_2D( 0, 0, 0, 0 ) 267 274 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( pssh(ji+1,jj ,Kmm) - pssh(ji ,jj ,Kmm) ) * r1_e1u(ji,jj) 268 275 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( pssh(ji ,jj+1,Kmm) - pssh(ji ,jj ,Kmm) ) * r1_e2v(ji,jj) … … 272 279 ENDIF 273 280 ! 274 DO_2D _00_00281 DO_2D( 0, 0, 0, 0 ) ! Remove coriolis term (and possibly spg) from barotropic trend 275 282 zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj) 276 283 zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj) … … 284 291 IF( ln_apr_dyn ) THEN 285 292 IF( ln_bt_fw ) THEN ! FORWARD integration: use kt+1/2 pressure (NOW+1/2) 286 DO_2D _00_00293 DO_2D( 0, 0, 0, 0 ) 287 294 zu_frc(ji,jj) = zu_frc(ji,jj) + grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj) 288 295 zv_frc(ji,jj) = zv_frc(ji,jj) + grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj) … … 290 297 ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 pressure (NOW) 291 298 zztmp = grav * r1_2 292 DO_2D _00_00299 DO_2D( 0, 0, 0, 0 ) 293 300 zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) & 294 301 & + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj) … … 302 309 ! ! ---------------------------------- ! 303 310 IF( ln_bt_fw ) THEN ! Add wind forcing 304 DO_2D _00_00311 DO_2D( 0, 0, 0, 0 ) 305 312 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rho0 * utau(ji,jj) * r1_hu(ji,jj,Kmm) 306 313 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rho0 * vtau(ji,jj) * r1_hv(ji,jj,Kmm) … … 308 315 ELSE 309 316 zztmp = r1_rho0 * r1_2 310 DO_2D _00_00317 DO_2D( 0, 0, 0, 0 ) 311 318 zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu(ji,jj,Kmm) 312 319 zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv(ji,jj,Kmm) … … 468 475 ! 469 476 ! ! ocean u- and v-depth at mid-step (separate DO-loops remove the need of a lbc_lnk) 470 DO_2D _11_10477 DO_2D( 1, 1, 1, 0 ) ! not jpi-column 471 478 zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj) & 472 479 & * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 473 480 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj) 474 481 END_2D 475 DO_2D _10_11482 DO_2D( 1, 0, 1, 1 ) ! not jpj-row 476 483 zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * r1_e1e2v(ji,jj) & 477 484 & * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & … … 508 515 !-- ssh = ssh - delta_t' * [ frc + div( flux ) ] --! 509 516 !-------------------------------------------------------------------------! 510 DO_2D _00_00517 DO_2D( 0, 0, 0, 0 ) 511 518 zhdiv = ( zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1) ) * r1_e1e2t(ji,jj) 512 519 ssha_e(ji,jj) = ( sshn_e(ji,jj) - rDt_e * ( zssh_frc(ji,jj) + zhdiv ) ) * ssmask(ji,jj) … … 514 521 ! 515 522 CALL lbc_lnk_multi( 'dynspg_ts', ssha_e, 'T', 1._wp, zhU, 'U', -1._wp, zhV, 'V', -1._wp ) 523 ! 524 ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T) 525 IF( ln_bdy ) CALL bdy_ssh( ssha_e ) 526 #if defined key_agrif 527 IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn ) 528 #endif 516 529 ! 517 530 ! ! Sum over sub-time-steps to compute advective velocities … … 525 538 END IF 526 539 ! 527 ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T)528 IF( ln_bdy ) CALL bdy_ssh( ssha_e )529 #if defined key_agrif530 IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )531 #endif532 540 ! 533 541 ! Sea Surface Height at u-,v-points (vvl case only) 534 542 IF( .NOT.ln_linssh ) THEN 535 DO_2D _00_00543 DO_2D( 0, 0, 0, 0 ) 536 544 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 537 545 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & … … 553 561 ! ! Surface pressure gradient 554 562 zldg = ( 1._wp - rn_scal_load ) * grav ! local factor 555 DO_2D _00_00563 DO_2D( 0, 0, 0, 0 ) 556 564 zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 557 565 zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) … … 567 575 ! at each time step. We however keep them constant here for optimization. 568 576 ! Recall that zhU and zhV hold fluxes at jn+0.5 (extrapolated not backward interpolated) 569 CALL dyn_cor_2d( zh up2_e, zhvp2_e, ua_e, va_e, zhU, zhV, zu_trd, zv_trd )577 CALL dyn_cor_2d( zhtp2_e, zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV, zu_trd, zv_trd ) 570 578 ! 571 579 ! Add tidal astronomical forcing if defined 572 580 IF ( ln_tide .AND. ln_tide_pot ) THEN 573 DO_2D _00_00581 DO_2D( 0, 0, 0, 0 ) 574 582 zu_trd(ji,jj) = zu_trd(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 575 583 zv_trd(ji,jj) = zv_trd(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) … … 580 588 !jth do implicitly instead 581 589 IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs 582 DO_2D _00_00590 DO_2D( 0, 0, 0, 0 ) 583 591 zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj) 584 592 zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj) … … 598 606 !------------------------------------------------------------------------------------------------------------------------! 599 607 IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form 600 DO_2D _00_00608 DO_2D( 0, 0, 0, 0 ) 601 609 ua_e(ji,jj) = ( un_e(ji,jj) & 602 610 & + rDt_e * ( zu_spg(ji,jj) & … … 613 621 ! 614 622 ELSE !* Flux form 615 DO_2D _00_00623 DO_2D( 0, 0, 0, 0 ) 616 624 ! ! hu_e, hv_e hold depth at jn, zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2 617 625 ! ! backward interpolated depth used in spg terms at jn+1/2 … … 637 645 !jth implicit bottom friction: 638 646 IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs 639 DO_2D _00_00647 DO_2D( 0, 0, 0, 0 ) 640 648 ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rDt_e * zCdU_u(ji,jj) * hur_e(ji,jj)) 641 649 va_e(ji,jj) = va_e(ji,jj) /(1.0 - rDt_e * zCdU_v(ji,jj) * hvr_e(ji,jj)) … … 643 651 ENDIF 644 652 645 IF( .NOT.ln_linssh ) THEN 653 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 646 654 hu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1) 647 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) ) 648 656 hv_e (2:jpim1,2:jpjm1) = hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1) 649 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) ) 658 ENDIF 659 ! 660 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 650 661 CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp & 651 662 & , hu_e , 'U', 1._wp, hv_e , 'V', 1._wp & … … 654 665 CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp ) 655 666 ENDIF 656 !657 !658 667 ! ! open boundaries 659 668 IF( ln_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e ) … … 703 712 IF (ln_bt_fw) THEN 704 713 IF( .NOT.( kt == nit000 .AND. l_1st_euler ) ) THEN 705 DO_2D _11_11714 DO_2D( 1, 1, 1, 1 ) 706 715 zun_save = un_adv(ji,jj) 707 716 zvn_save = vn_adv(ji,jj) … … 734 743 ELSE 735 744 ! At this stage, pssh(:,:,:,Krhs) has been corrected: compute new depths at velocity points 736 DO_2D _10_10745 DO_2D( 1, 0, 1, 0 ) 737 746 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 738 747 & * ( e1e2t(ji ,jj) * pssh(ji ,jj,Kaa) & … … 891 900 ! ! --------------- 892 901 IF( ln_rstart .AND. ln_bt_fw .AND. (.NOT.l_1st_euler) ) THEN !* Read the restart file 893 CALL iom_get( numror, jpdom_auto glo, 'ub2_b' , ub2_b (:,:), ldxios = lrxios )894 CALL iom_get( numror, jpdom_auto glo, 'vb2_b' , vb2_b (:,:), ldxios = lrxios )895 CALL iom_get( numror, jpdom_auto glo, 'un_bf' , un_bf (:,:), ldxios = lrxios )896 CALL iom_get( numror, jpdom_auto glo, 'vn_bf' , vn_bf (:,:), ldxios = lrxios )902 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 ) 897 906 IF( .NOT.ln_bt_av ) THEN 898 CALL iom_get( numror, jpdom_auto glo, 'sshbb_e' , sshbb_e(:,:), ldxios = lrxios )899 CALL iom_get( numror, jpdom_auto glo, 'ubb_e' , ubb_e(:,:), ldxios = lrxios )900 CALL iom_get( numror, jpdom_auto glo, 'vbb_e' , vbb_e(:,:), ldxios = lrxios )901 CALL iom_get( numror, jpdom_auto glo, 'sshb_e' , sshb_e(:,:), ldxios = lrxios )902 CALL iom_get( numror, jpdom_auto glo, 'ub_e' , ub_e(:,:), ldxios = lrxios )903 CALL iom_get( numror, jpdom_auto glo, 'vb_e' , vb_e(:,:), ldxios = lrxios )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 ) 904 913 ENDIF 905 914 #if defined key_agrif 906 915 ! Read time integrated fluxes 907 916 IF ( .NOT.Agrif_Root() ) THEN 908 CALL iom_get( numror, jpdom_auto glo, 'ub2_i_b' , ub2_i_b(:,:), ldxios = lrxios )909 CALL iom_get( numror, jpdom_auto glo, 'vb2_i_b' , vb2_i_b(:,:), ldxios = lrxios )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 ) 910 919 ENDIF 911 920 #endif … … 966 975 ! Max courant number for ext. grav. waves 967 976 ! 968 DO_2D _11_11977 DO_2D( 0, 0, 0, 0 ) 969 978 zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 970 979 zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj) … … 972 981 END_2D 973 982 ! 974 zcmax = MAXVAL( zcu( :,:) )983 zcmax = MAXVAL( zcu(Nis0:Nie0,Njs0:Nje0) ) 975 984 CALL mpp_max( 'dynspg_ts', zcmax ) 976 985 … … 1088 1097 ! 1089 1098 SELECT CASE( nvor_scheme ) 1090 CASE( np_EEN ) != EEN scheme using e3f (energy & enstrophy scheme)1099 CASE( np_EEN ) != EEN scheme using e3f energy & enstrophy scheme 1091 1100 SELECT CASE( nn_een_e3f ) !* ff_f/e3 at F-point 1092 1101 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 1093 DO_2D _10_101102 DO_2D( 1, 0, 1, 0 ) 1094 1103 zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + & 1095 1104 & ht(ji ,jj ) + ht(ji+1,jj ) ) * 0.25_wp … … 1097 1106 END_2D 1098 1107 CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask) 1099 DO_2D _10_101108 DO_2D( 1, 0, 1, 0 ) 1100 1109 zwz(ji,jj) = ( ht (ji ,jj+1) + ht (ji+1,jj+1) & 1101 1110 & + ht (ji ,jj ) + ht (ji+1,jj ) ) & … … 1108 1117 ! 1109 1118 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 1110 DO_2D _01_011119 DO_2D( 0, 1, 0, 1 ) 1111 1120 ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 1112 1121 ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) … … 1115 1124 END_2D 1116 1125 ! 1117 CASE( np_EET ) != EEN scheme using e3t (energy conserving scheme)1126 CASE( np_EET ) != EEN scheme using e3t energy conserving scheme 1118 1127 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 1119 DO_2D _01_011128 DO_2D( 0, 1, 0, 1 ) 1120 1129 z1_ht = ssmask(ji,jj) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) ) 1121 1130 ftne(ji,jj) = ( ff_f(ji-1,jj ) + ff_f(ji ,jj ) + ff_f(ji ,jj-1) ) * z1_ht … … 1150 1159 ! 1151 1160 !zhf(:,:) = hbatf(:,:) 1152 DO_2D _10_101161 DO_2D( 1, 0, 1, 0 ) 1153 1162 zhf(ji,jj) = ( ht_0 (ji,jj ) + ht_0 (ji+1,jj ) & 1154 1163 & + ht_0 (ji,jj+1) + ht_0 (ji+1,jj+1) ) & … … 1169 1178 CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp ) 1170 1179 ! JC: TBC. hf should be greater than 0 1171 DO_2D _11_111180 DO_2D( 1, 1, 1, 1 ) 1172 1181 IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj) 1173 1182 END_2D … … 1179 1188 1180 1189 1181 SUBROUTINE dyn_cor_2d( ph u, phv, punb, pvnb, zhU, zhV, zu_trd, zv_trd )1190 SUBROUTINE dyn_cor_2d( pht, phu, phv, punb, pvnb, zhU, zhV, zu_trd, zv_trd ) 1182 1191 !!--------------------------------------------------------------------- 1183 1192 !! *** ROUTINE dyn_cor_2d *** … … 1187 1196 INTEGER :: ji ,jj ! dummy loop indices 1188 1197 REAL(wp) :: zx1, zx2, zy1, zy2, z1_hu, z1_hv ! - - 1189 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: ph u, phv, punb, pvnb, zhU, zhV1198 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pht, phu, phv, punb, pvnb, zhU, zhV 1190 1199 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: zu_trd, zv_trd 1191 1200 !!---------------------------------------------------------------------- 1192 1201 SELECT CASE( nvor_scheme ) 1193 1202 CASE( np_ENT ) ! enstrophy conserving scheme (f-point) 1194 DO_2D _00_001203 DO_2D( 0, 0, 0, 0 ) 1195 1204 z1_hu = ssumask(ji,jj) / ( phu(ji,jj) + 1._wp - ssumask(ji,jj) ) 1196 1205 z1_hv = ssvmask(ji,jj) / ( phv(ji,jj) + 1._wp - ssvmask(ji,jj) ) 1197 1206 zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu & 1198 & * ( e1e2t(ji+1,jj)* ht(ji+1,jj)*ff_t(ji+1,jj) * ( pvnb(ji+1,jj) + pvnb(ji+1,jj-1) ) &1199 & + e1e2t(ji ,jj)* ht(ji ,jj)*ff_t(ji ,jj) * ( pvnb(ji ,jj) + pvnb(ji ,jj-1) ) )1207 & * ( e1e2t(ji+1,jj)*pht(ji+1,jj)*ff_t(ji+1,jj) * ( pvnb(ji+1,jj) + pvnb(ji+1,jj-1) ) & 1208 & + e1e2t(ji ,jj)*pht(ji ,jj)*ff_t(ji ,jj) * ( pvnb(ji ,jj) + pvnb(ji ,jj-1) ) ) 1200 1209 ! 1201 1210 zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv & 1202 & * ( e1e2t(ji,jj+1)* ht(ji,jj+1)*ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) ) &1203 & + e1e2t(ji,jj )* ht(ji,jj )*ff_t(ji,jj ) * ( punb(ji,jj ) + punb(ji-1,jj ) ) )1211 & * ( e1e2t(ji,jj+1)*pht(ji,jj+1)*ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) ) & 1212 & + e1e2t(ji,jj )*pht(ji,jj )*ff_t(ji,jj ) * ( punb(ji,jj ) + punb(ji-1,jj ) ) ) 1204 1213 END_2D 1205 1214 ! 1206 1215 CASE( np_ENE , np_MIX ) ! energy conserving scheme (t-point) ENE or MIX 1207 DO_2D _00_001216 DO_2D( 0, 0, 0, 0 ) 1208 1217 zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj) 1209 1218 zy2 = ( zhV(ji,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj) … … 1216 1225 ! 1217 1226 CASE( np_ENS ) ! enstrophy conserving scheme (f-point) 1218 DO_2D _00_001227 DO_2D( 0, 0, 0, 0 ) 1219 1228 zy1 = r1_8 * ( zhV(ji ,jj-1) + zhV(ji+1,jj-1) & 1220 1229 & + zhV(ji ,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj) … … 1226 1235 ! 1227 1236 CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f) 1228 DO_2D _00_001237 DO_2D( 0, 0, 0, 0 ) 1229 1238 zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zhV(ji ,jj ) & 1230 1239 & + ftnw(ji+1,jj) * zhV(ji+1,jj ) & … … 1260 1269 ! 1261 1270 IF( ln_wd_dl_rmp ) THEN 1262 DO_2D _11_111271 DO_2D( 1, 1, 1, 1 ) 1263 1272 IF ( pssh(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN 1264 1273 ! IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN … … 1271 1280 END_2D 1272 1281 ELSE 1273 DO_2D _11_111282 DO_2D( 1, 1, 1, 1 ) 1274 1283 IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN ; ptmsk(ji,jj) = 1._wp 1275 1284 ELSE ; ptmsk(ji,jj) = 0._wp … … 1299 1308 !!---------------------------------------------------------------------- 1300 1309 ! 1301 DO_2D _11_101310 DO_2D( 1, 1, 1, 0 ) ! not jpi-column 1302 1311 IF ( phU(ji,jj) > 0._wp ) THEN ; pUmsk(ji,jj) = pTmsk(ji ,jj) 1303 1312 ELSE ; pUmsk(ji,jj) = pTmsk(ji+1,jj) … … 1307 1316 END_2D 1308 1317 ! 1309 DO_2D _10_111318 DO_2D( 1, 0, 1, 1 ) ! not jpj-row 1310 1319 IF ( phV(ji,jj) > 0._wp ) THEN ; pVmsk(ji,jj) = pTmsk(ji,jj ) 1311 1320 ELSE ; pVmsk(ji,jj) = pTmsk(ji,jj+1) … … 1329 1338 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy 1330 1339 !!---------------------------------------------------------------------- 1331 DO_2D _00_001340 DO_2D( 0, 0, 0, 0 ) 1332 1341 ll_tmp1 = MIN( pshn(ji,jj) , pshn(ji+1,jj) ) > & 1333 1342 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & … … 1396 1405 ! !== Set the barotropic drag coef. ==! 1397 1406 ! 1398 IF( ln_isfcav ) THEN ! top+bottom friction (ocean cavities)1407 IF( ln_isfcav.OR.ln_drgice_imp ) THEN ! top+bottom friction (ocean cavities) 1399 1408 1400 DO_2D _00_001409 DO_2D( 0, 0, 0, 0 ) 1401 1410 pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) 1402 1411 pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) 1403 1412 END_2D 1404 1413 ELSE ! bottom friction only 1405 DO_2D _00_001414 DO_2D( 0, 0, 0, 0 ) 1406 1415 pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) 1407 1416 pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) … … 1413 1422 IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW bottom baroclinic velocities 1414 1423 1415 DO_2D _00_001424 DO_2D( 0, 0, 0, 0 ) 1416 1425 ikbu = mbku(ji,jj) 1417 1426 ikbv = mbkv(ji,jj) … … 1421 1430 ELSE ! CENTRED integration: use BEFORE bottom baroclinic velocities 1422 1431 1423 DO_2D _00_001432 DO_2D( 0, 0, 0, 0 ) 1424 1433 ikbu = mbku(ji,jj) 1425 1434 ikbv = mbkv(ji,jj) … … 1431 1440 IF( ln_wd_il ) THEN ! W/D : use the "clipped" bottom friction !!gm explain WHY, please ! 1432 1441 zztmp = -1._wp / rDt_e 1433 DO_2D _00_001442 DO_2D( 0, 0, 0, 0 ) 1434 1443 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( & 1435 1444 & r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp ) … … 1439 1448 ELSE ! use "unclipped" drag (even if explicit friction is used in 3D calculation) 1440 1449 1441 DO_2D _00_001450 DO_2D( 0, 0, 0, 0 ) 1442 1451 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zu_i(ji,jj) 1443 1452 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zv_i(ji,jj) … … 1447 1456 ! !== TOP stress contribution from baroclinic velocities ==! (no W/D case) 1448 1457 ! 1449 IF( ln_isfcav ) THEN1458 IF( ln_isfcav.OR.ln_drgice_imp ) THEN 1450 1459 ! 1451 1460 IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW top baroclinic velocity 1452 1461 1453 DO_2D _00_001462 DO_2D( 0, 0, 0, 0 ) 1454 1463 iktu = miku(ji,jj) 1455 1464 iktv = mikv(ji,jj) … … 1459 1468 ELSE ! CENTRED integration: use BEFORE top baroclinic velocity 1460 1469 1461 DO_2D _00_001470 DO_2D( 0, 0, 0, 0 ) 1462 1471 iktu = miku(ji,jj) 1463 1472 iktv = mikv(ji,jj) … … 1469 1478 ! ! use "unclipped" top drag (even if explicit friction is used in 3D calculation) 1470 1479 1471 DO_2D _00_001480 DO_2D( 0, 0, 0, 0 ) 1472 1481 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zu_i(ji,jj) 1473 1482 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zv_i(ji,jj)
Note: See TracChangeset
for help on using the changeset viewer.
