- Timestamp:
- 2020-11-27T17:26:33+01:00 (4 years ago)
- Location:
- NEMO/branches/2020/tickets_icb_1900
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/tickets_icb_1900
- Property svn:externals
-
NEMO/branches/2020/tickets_icb_1900/src/OCE/DYN/dynspg_ts.F90
r13237 r13899 264 264 IF( ln_wd_il ) THEN ! W/D : limiter applied to spgspg 265 265 CALL wad_spg( pssh(:,:,Kmm), zcpx, zcpy ) ! Calculating W/D gravity filters, zcpx and zcpy 266 DO_2D _00_00266 DO_2D( 0, 0, 0, 0 ) ! SPG with the application of W/D gravity filters 267 267 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( pssh(ji+1,jj ,Kmm) - pssh(ji ,jj ,Kmm) ) & 268 268 & * r1_e1u(ji,jj) * zcpx(ji,jj) * wdrampu(ji,jj) !jth … … 271 271 END_2D 272 272 ELSE ! now suface pressure gradient 273 DO_2D _00_00273 DO_2D( 0, 0, 0, 0 ) 274 274 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( pssh(ji+1,jj ,Kmm) - pssh(ji ,jj ,Kmm) ) * r1_e1u(ji,jj) 275 275 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( pssh(ji ,jj+1,Kmm) - pssh(ji ,jj ,Kmm) ) * r1_e2v(ji,jj) … … 279 279 ENDIF 280 280 ! 281 DO_2D _00_00281 DO_2D( 0, 0, 0, 0 ) ! Remove coriolis term (and possibly spg) from barotropic trend 282 282 zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj) 283 283 zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj) … … 291 291 IF( ln_apr_dyn ) THEN 292 292 IF( ln_bt_fw ) THEN ! FORWARD integration: use kt+1/2 pressure (NOW+1/2) 293 DO_2D _00_00293 DO_2D( 0, 0, 0, 0 ) 294 294 zu_frc(ji,jj) = zu_frc(ji,jj) + grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj) 295 295 zv_frc(ji,jj) = zv_frc(ji,jj) + grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj) … … 297 297 ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 pressure (NOW) 298 298 zztmp = grav * r1_2 299 DO_2D _00_00299 DO_2D( 0, 0, 0, 0 ) 300 300 zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) & 301 301 & + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj) … … 309 309 ! ! ---------------------------------- ! 310 310 IF( ln_bt_fw ) THEN ! Add wind forcing 311 DO_2D _00_00311 DO_2D( 0, 0, 0, 0 ) 312 312 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rho0 * utau(ji,jj) * r1_hu(ji,jj,Kmm) 313 313 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rho0 * vtau(ji,jj) * r1_hv(ji,jj,Kmm) … … 315 315 ELSE 316 316 zztmp = r1_rho0 * r1_2 317 DO_2D _00_00317 DO_2D( 0, 0, 0, 0 ) 318 318 zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu(ji,jj,Kmm) 319 319 zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv(ji,jj,Kmm) … … 475 475 ! 476 476 ! ! ocean u- and v-depth at mid-step (separate DO-loops remove the need of a lbc_lnk) 477 DO_2D _11_10477 DO_2D( 1, 1, 1, 0 ) ! not jpi-column 478 478 zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj) & 479 479 & * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 480 480 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj) 481 481 END_2D 482 DO_2D _10_11482 DO_2D( 1, 0, 1, 1 ) ! not jpj-row 483 483 zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * r1_e1e2v(ji,jj) & 484 484 & * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & … … 515 515 !-- ssh = ssh - delta_t' * [ frc + div( flux ) ] --! 516 516 !-------------------------------------------------------------------------! 517 DO_2D _00_00517 DO_2D( 0, 0, 0, 0 ) 518 518 zhdiv = ( zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1) ) * r1_e1e2t(ji,jj) 519 519 ssha_e(ji,jj) = ( sshn_e(ji,jj) - rDt_e * ( zssh_frc(ji,jj) + zhdiv ) ) * ssmask(ji,jj) 520 520 END_2D 521 ! 522 CALL lbc_lnk_multi( 'dynspg_ts', ssha_e, 'T', 1._wp, zhU, 'U', -1._wp, zhV, 'V', -1._wp ) 521 523 ! 522 524 ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T) … … 525 527 IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn ) 526 528 #endif 527 528 CALL lbc_lnk_multi( 'dynspg_ts', ssha_e, 'T', 1._wp, zhU, 'U', -1._wp, zhV, 'V', -1._wp )529 529 ! 530 530 ! ! Sum over sub-time-steps to compute advective velocities … … 541 541 ! Sea Surface Height at u-,v-points (vvl case only) 542 542 IF( .NOT.ln_linssh ) THEN 543 DO_2D _00_00543 DO_2D( 0, 0, 0, 0 ) 544 544 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 545 545 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & … … 561 561 ! ! Surface pressure gradient 562 562 zldg = ( 1._wp - rn_scal_load ) * grav ! local factor 563 DO_2D _00_00563 DO_2D( 0, 0, 0, 0 ) 564 564 zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 565 565 zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) … … 579 579 ! Add tidal astronomical forcing if defined 580 580 IF ( ln_tide .AND. ln_tide_pot ) THEN 581 DO_2D _00_00581 DO_2D( 0, 0, 0, 0 ) 582 582 zu_trd(ji,jj) = zu_trd(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 583 583 zv_trd(ji,jj) = zv_trd(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) … … 588 588 !jth do implicitly instead 589 589 IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs 590 DO_2D _00_00590 DO_2D( 0, 0, 0, 0 ) 591 591 zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj) 592 592 zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj) … … 606 606 !------------------------------------------------------------------------------------------------------------------------! 607 607 IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form 608 DO_2D _00_00608 DO_2D( 0, 0, 0, 0 ) 609 609 ua_e(ji,jj) = ( un_e(ji,jj) & 610 610 & + rDt_e * ( zu_spg(ji,jj) & … … 621 621 ! 622 622 ELSE !* Flux form 623 DO_2D _00_00623 DO_2D( 0, 0, 0, 0 ) 624 624 ! ! hu_e, hv_e hold depth at jn, zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2 625 625 ! ! backward interpolated depth used in spg terms at jn+1/2 … … 645 645 !jth implicit bottom friction: 646 646 IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs 647 DO_2D _00_00647 DO_2D( 0, 0, 0, 0 ) 648 648 ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rDt_e * zCdU_u(ji,jj) * hur_e(ji,jj)) 649 649 va_e(ji,jj) = va_e(ji,jj) /(1.0 - rDt_e * zCdU_v(ji,jj) * hvr_e(ji,jj)) … … 657 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 658 ENDIF 659 ! ! open boundaries660 IF( ln_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )661 #if defined key_agrif662 IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jn ) ! Agrif663 #endif664 659 ! 665 660 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) … … 670 665 CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp ) 671 666 ENDIF 672 ! 667 ! ! open boundaries 668 IF( ln_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e ) 669 #if defined key_agrif 670 IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jn ) ! Agrif 671 #endif 673 672 ! !* Swap 674 673 ! ! ---- … … 713 712 IF (ln_bt_fw) THEN 714 713 IF( .NOT.( kt == nit000 .AND. l_1st_euler ) ) THEN 715 DO_2D _11_11714 DO_2D( 1, 1, 1, 1 ) 716 715 zun_save = un_adv(ji,jj) 717 716 zvn_save = vn_adv(ji,jj) … … 744 743 ELSE 745 744 ! At this stage, pssh(:,:,:,Krhs) has been corrected: compute new depths at velocity points 746 DO_2D _10_10745 DO_2D( 1, 0, 1, 0 ) 747 746 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 748 747 & * ( e1e2t(ji ,jj) * pssh(ji ,jj,Kaa) & … … 901 900 ! ! --------------- 902 901 IF( ln_rstart .AND. ln_bt_fw .AND. (.NOT.l_1st_euler) ) THEN !* Read the restart file 903 CALL iom_get( numror, jpdom_auto glo, 'ub2_b' , ub2_b (:,:), ldxios = lrxios )904 CALL iom_get( numror, jpdom_auto glo, 'vb2_b' , vb2_b (:,:), ldxios = lrxios )905 CALL iom_get( numror, jpdom_auto glo, 'un_bf' , un_bf (:,:), ldxios = lrxios )906 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 ) 907 906 IF( .NOT.ln_bt_av ) THEN 908 CALL iom_get( numror, jpdom_auto glo, 'sshbb_e' , sshbb_e(:,:), ldxios = lrxios )909 CALL iom_get( numror, jpdom_auto glo, 'ubb_e' , ubb_e(:,:), ldxios = lrxios )910 CALL iom_get( numror, jpdom_auto glo, 'vbb_e' , vbb_e(:,:), ldxios = lrxios )911 CALL iom_get( numror, jpdom_auto glo, 'sshb_e' , sshb_e(:,:), ldxios = lrxios )912 CALL iom_get( numror, jpdom_auto glo, 'ub_e' , ub_e(:,:), ldxios = lrxios )913 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 ) 914 913 ENDIF 915 914 #if defined key_agrif 916 915 ! Read time integrated fluxes 917 916 IF ( .NOT.Agrif_Root() ) THEN 918 CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b' , ub2_i_b(:,:), ldxios = lrxios ) 919 CALL iom_get( numror, jpdom_autoglo, '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 ) 919 ELSE 920 ub2_i_b(:,:) = 0._wp ; vb2_i_b(:,:) = 0._wp ! used in the 1st update of agrif 920 921 ENDIF 921 922 #endif … … 923 924 IF(lwp) WRITE(numout,*) 924 925 IF(lwp) WRITE(numout,*) ' ==>>> start from rest: set barotropic values to 0' 925 ub2_b (:,:) = 0._wp ; vb2_b(:,:) = 0._wp ! used in the 1st interpol of agrif926 un_adv (:,:) = 0._wp ; vn_adv(:,:) = 0._wp ! used in the 1st interpol of agrif927 un_bf (:,:) = 0._wp ; vn_bf(:,:) = 0._wp ! used in the 1st update of agrif926 ub2_b (:,:) = 0._wp ; vb2_b (:,:) = 0._wp ! used in the 1st interpol of agrif 927 un_adv (:,:) = 0._wp ; vn_adv (:,:) = 0._wp ! used in the 1st interpol of agrif 928 un_bf (:,:) = 0._wp ; vn_bf (:,:) = 0._wp ! used in the 1st update of agrif 928 929 #if defined key_agrif 929 IF ( .NOT.Agrif_Root() ) THEN 930 ub2_i_b(:,:) = 0._wp ; vb2_i_b(:,:) = 0._wp ! used in the 1st update of agrif 931 ENDIF 930 ub2_i_b(:,:) = 0._wp ; vb2_i_b(:,:) = 0._wp ! used in the 1st update of agrif 932 931 #endif 933 932 ENDIF … … 976 975 ! Max courant number for ext. grav. waves 977 976 ! 978 DO_2D _11_11977 DO_2D( 0, 0, 0, 0 ) 979 978 zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 980 979 zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj) … … 982 981 END_2D 983 982 ! 984 zcmax = MAXVAL( zcu( :,:) )983 zcmax = MAXVAL( zcu(Nis0:Nie0,Njs0:Nje0) ) 985 984 CALL mpp_max( 'dynspg_ts', zcmax ) 986 985 … … 1101 1100 SELECT CASE( nn_een_e3f ) !* ff_f/e3 at F-point 1102 1101 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 1103 DO_2D _10_101102 DO_2D( 1, 0, 1, 0 ) 1104 1103 zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + & 1105 1104 & ht(ji ,jj ) + ht(ji+1,jj ) ) * 0.25_wp … … 1107 1106 END_2D 1108 1107 CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask) 1109 DO_2D _10_101108 DO_2D( 1, 0, 1, 0 ) 1110 1109 zwz(ji,jj) = ( ht (ji ,jj+1) + ht (ji+1,jj+1) & 1111 1110 & + ht (ji ,jj ) + ht (ji+1,jj ) ) & … … 1118 1117 ! 1119 1118 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 1120 DO_2D _01_011119 DO_2D( 0, 1, 0, 1 ) 1121 1120 ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 1122 1121 ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) … … 1127 1126 CASE( np_EET ) != EEN scheme using e3t energy conserving scheme 1128 1127 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 1129 DO_2D _01_011128 DO_2D( 0, 1, 0, 1 ) 1130 1129 z1_ht = ssmask(ji,jj) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) ) 1131 1130 ftne(ji,jj) = ( ff_f(ji-1,jj ) + ff_f(ji ,jj ) + ff_f(ji ,jj-1) ) * z1_ht … … 1160 1159 ! 1161 1160 !zhf(:,:) = hbatf(:,:) 1162 DO_2D _10_101161 DO_2D( 1, 0, 1, 0 ) 1163 1162 zhf(ji,jj) = ( ht_0 (ji,jj ) + ht_0 (ji+1,jj ) & 1164 1163 & + ht_0 (ji,jj+1) + ht_0 (ji+1,jj+1) ) & … … 1179 1178 CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp ) 1180 1179 ! JC: TBC. hf should be greater than 0 1181 DO_2D _11_111180 DO_2D( 1, 1, 1, 1 ) 1182 1181 IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj) 1183 1182 END_2D … … 1202 1201 SELECT CASE( nvor_scheme ) 1203 1202 CASE( np_ENT ) ! enstrophy conserving scheme (f-point) 1204 DO_2D _00_001203 DO_2D( 0, 0, 0, 0 ) 1205 1204 z1_hu = ssumask(ji,jj) / ( phu(ji,jj) + 1._wp - ssumask(ji,jj) ) 1206 1205 z1_hv = ssvmask(ji,jj) / ( phv(ji,jj) + 1._wp - ssvmask(ji,jj) ) … … 1215 1214 ! 1216 1215 CASE( np_ENE , np_MIX ) ! energy conserving scheme (t-point) ENE or MIX 1217 DO_2D _00_001216 DO_2D( 0, 0, 0, 0 ) 1218 1217 zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj) 1219 1218 zy2 = ( zhV(ji,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj) … … 1226 1225 ! 1227 1226 CASE( np_ENS ) ! enstrophy conserving scheme (f-point) 1228 DO_2D _00_001227 DO_2D( 0, 0, 0, 0 ) 1229 1228 zy1 = r1_8 * ( zhV(ji ,jj-1) + zhV(ji+1,jj-1) & 1230 1229 & + zhV(ji ,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj) … … 1236 1235 ! 1237 1236 CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f) 1238 DO_2D _00_001237 DO_2D( 0, 0, 0, 0 ) 1239 1238 zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zhV(ji ,jj ) & 1240 1239 & + ftnw(ji+1,jj) * zhV(ji+1,jj ) & … … 1270 1269 ! 1271 1270 IF( ln_wd_dl_rmp ) THEN 1272 DO_2D _11_111271 DO_2D( 1, 1, 1, 1 ) 1273 1272 IF ( pssh(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN 1274 1273 ! IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN … … 1281 1280 END_2D 1282 1281 ELSE 1283 DO_2D _11_111282 DO_2D( 1, 1, 1, 1 ) 1284 1283 IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN ; ptmsk(ji,jj) = 1._wp 1285 1284 ELSE ; ptmsk(ji,jj) = 0._wp … … 1309 1308 !!---------------------------------------------------------------------- 1310 1309 ! 1311 DO_2D _11_101310 DO_2D( 1, 1, 1, 0 ) ! not jpi-column 1312 1311 IF ( phU(ji,jj) > 0._wp ) THEN ; pUmsk(ji,jj) = pTmsk(ji ,jj) 1313 1312 ELSE ; pUmsk(ji,jj) = pTmsk(ji+1,jj) … … 1317 1316 END_2D 1318 1317 ! 1319 DO_2D _10_111318 DO_2D( 1, 0, 1, 1 ) ! not jpj-row 1320 1319 IF ( phV(ji,jj) > 0._wp ) THEN ; pVmsk(ji,jj) = pTmsk(ji,jj ) 1321 1320 ELSE ; pVmsk(ji,jj) = pTmsk(ji,jj+1) … … 1339 1338 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy 1340 1339 !!---------------------------------------------------------------------- 1341 DO_2D _00_001340 DO_2D( 0, 0, 0, 0 ) 1342 1341 ll_tmp1 = MIN( pshn(ji,jj) , pshn(ji+1,jj) ) > & 1343 1342 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & … … 1406 1405 ! !== Set the barotropic drag coef. ==! 1407 1406 ! 1408 IF( ln_isfcav ) THEN ! top+bottom friction (ocean cavities)1407 IF( ln_isfcav.OR.ln_drgice_imp ) THEN ! top+bottom friction (ocean cavities) 1409 1408 1410 DO_2D _00_001409 DO_2D( 0, 0, 0, 0 ) 1411 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) ) 1412 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) ) 1413 1412 END_2D 1414 1413 ELSE ! bottom friction only 1415 DO_2D _00_001414 DO_2D( 0, 0, 0, 0 ) 1416 1415 pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) 1417 1416 pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) … … 1423 1422 IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW bottom baroclinic velocities 1424 1423 1425 DO_2D _00_001424 DO_2D( 0, 0, 0, 0 ) 1426 1425 ikbu = mbku(ji,jj) 1427 1426 ikbv = mbkv(ji,jj) … … 1431 1430 ELSE ! CENTRED integration: use BEFORE bottom baroclinic velocities 1432 1431 1433 DO_2D _00_001432 DO_2D( 0, 0, 0, 0 ) 1434 1433 ikbu = mbku(ji,jj) 1435 1434 ikbv = mbkv(ji,jj) … … 1441 1440 IF( ln_wd_il ) THEN ! W/D : use the "clipped" bottom friction !!gm explain WHY, please ! 1442 1441 zztmp = -1._wp / rDt_e 1443 DO_2D _00_001442 DO_2D( 0, 0, 0, 0 ) 1444 1443 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( & 1445 1444 & r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp ) … … 1449 1448 ELSE ! use "unclipped" drag (even if explicit friction is used in 3D calculation) 1450 1449 1451 DO_2D _00_001450 DO_2D( 0, 0, 0, 0 ) 1452 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) 1453 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) … … 1457 1456 ! !== TOP stress contribution from baroclinic velocities ==! (no W/D case) 1458 1457 ! 1459 IF( ln_isfcav ) THEN1458 IF( ln_isfcav.OR.ln_drgice_imp ) THEN 1460 1459 ! 1461 1460 IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW top baroclinic velocity 1462 1461 1463 DO_2D _00_001462 DO_2D( 0, 0, 0, 0 ) 1464 1463 iktu = miku(ji,jj) 1465 1464 iktv = mikv(ji,jj) … … 1469 1468 ELSE ! CENTRED integration: use BEFORE top baroclinic velocity 1470 1469 1471 DO_2D _00_001470 DO_2D( 0, 0, 0, 0 ) 1472 1471 iktu = miku(ji,jj) 1473 1472 iktv = mikv(ji,jj) … … 1479 1478 ! ! use "unclipped" top drag (even if explicit friction is used in 3D calculation) 1480 1479 1481 DO_2D _00_001480 DO_2D( 0, 0, 0, 0 ) 1482 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) 1483 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.