Changeset 11234
- Timestamp:
- 2019-07-09T18:09:16+02:00 (6 years ago)
- Location:
- NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE
- Files:
-
- 3 edited
-
BDY/bdydyn3d.F90 (modified) (1 diff)
-
BDY/bdyvol.F90 (modified) (6 diffs)
-
DYN/dynspg_ts.F90 (modified) (27 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdydyn3d.F90
r11210 r11234 186 186 igrd = 3 ! Copying tangential velocity into bdy points 187 187 IF( llrim0 ) THEN ; ibeg = 1 ; iend = idx%nblenrim0(igrd) 188 ELSE ; ibeg = idx%nblenrim0(igrd)+1 ; iend = idx%nblenrim(igrd)188 ELSE ; ibeg = idx%nblenrim0(igrd)+1 ; iend = idx%nblenrim(igrd) 189 189 ENDIF 190 190 DO jb = ibeg, iend -
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdyvol.F90
r11049 r11234 99 99 ii = idx%nbi(jb,jgrd) 100 100 ij = idx%nbj(jb,jgrd) 101 IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE ! to remove101 IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE ! sum : else halo couted twice 102 102 zubtpecor = zubtpecor + idx%flagu(jb,jgrd) * pua2d(ii,ij) * e2u(ii,ij) * phu(ii,ij) * tmask_i(ii,ij) * tmask_i(ii+1,ij) 103 103 END DO … … 106 106 ii = idx%nbi(jb,jgrd) 107 107 ij = idx%nbj(jb,jgrd) 108 IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE ! to remove108 IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE ! sum : else halo couted twice 109 109 zubtpecor = zubtpecor + idx%flagv(jb,jgrd) * pva2d(ii,ij) * e1v(ii,ij) * phv(ii,ij) * tmask_i(ii,ij) * tmask_i(ii,ij+1) 110 110 END DO … … 128 128 ii = idx%nbi(jb,jgrd) 129 129 ij = idx%nbj(jb,jgrd) 130 IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE ! sum : else halo couted twice130 !IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE ! to remove ? 131 131 pua2d(ii,ij) = pua2d(ii,ij) - idx%flagu(jb,jgrd) * zubtpecor * tmask_i(ii,ij) * tmask_i(ii+1,ij) 132 132 END DO … … 135 135 ii = idx%nbi(jb,jgrd) 136 136 ij = idx%nbj(jb,jgrd) 137 IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE ! sum : else halo couted twice137 !IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE ! to remove ? 138 138 pva2d(ii,ij) = pva2d(ii,ij) - idx%flagv(jb,jgrd) * zubtpecor * tmask_i(ii,ij) * tmask_i(ii,ij+1) 139 139 END DO … … 154 154 ii = idx%nbi(jb,jgrd) 155 155 ij = idx%nbj(jb,jgrd) 156 IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE ! to remove?156 IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE 157 157 ztranst = ztranst + idx%flagu(jb,jgrd) * pua2d(ii,ij) * e2u(ii,ij) * phu(ii,ij) * tmask_i(ii,ij) * tmask_i(ii+1,ij) 158 158 END DO … … 161 161 ii = idx%nbi(jb,jgrd) 162 162 ij = idx%nbj(jb,jgrd) 163 IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE ! to remove?163 IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj ) CYCLE 164 164 ztranst = ztranst + idx%flagv(jb,jgrd) * pva2d(ii,ij) * e1v(ii,ij) * phv(ii,ij) * tmask_i(ii,ij) * tmask_i(ii,ij+1) 165 165 END DO -
NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DYN/dynspg_ts.F90
r11223 r11234 64 64 USE diatmb ! Top,middle,bottom output 65 65 66 USE iom ! to remove 67 66 68 IMPLICIT NONE 67 69 PRIVATE … … 104 106 ! 105 107 ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT=ierr(1) ) 106 !107 108 IF( ln_dynvor_een .OR. ln_dynvor_eeT ) & 108 & ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 109 & ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(2) ) 109 & ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , ftsw(jpi,jpj) , ftse(jpi,jpj), STAT=ierr(2) ) 110 110 ! 111 111 ALLOCATE( un_adv(jpi,jpj), vn_adv(jpi,jpj) , STAT=ierr(3) ) … … 153 153 INTEGER :: ikbv, iktv ! - - 154 154 REAL(wp) :: r1_2dt_b, z2dt_bf ! local scalars 155 REAL(wp) :: zx1, zx2, z u_spg, zhura, z1_hu ! - -156 REAL(wp) :: zy1, zy2, z v_spg, zhvra, z1_hv ! - -155 REAL(wp) :: zx1, zx2, zhura , z1_hu ! - - 156 REAL(wp) :: zy1, zy2, zhvra , z1_hv ! - - 157 157 REAL(wp) :: za0, za1, za2, za3 ! - - 158 REAL(wp) :: zmdi, zztmp , z1_ht ! - - 159 REAL(wp), DIMENSION(jpi,jpj) :: zsshp2_e, zhf 160 REAL(wp), DIMENSION(jpi,jpj) :: zwx, zu_trd, zu_frc, zssh_frc 161 REAL(wp), DIMENSION(jpi,jpj) :: zwy, zv_trd, zv_frc, zhdiv 158 REAL(wp) :: zmdi, zztmp, zldg , z1_ht ! - - 159 REAL(wp) :: zhu_bck, zhv_bck ! - - 160 REAL(wp) :: zun_save, zvn_save ! - - 161 REAL(wp), DIMENSION(jpi,jpj) :: zsshp2_e 162 REAL(wp), DIMENSION(jpi,jpj) :: zwx, zu_trd, zu_frc, zu_spg, zssh_frc 163 REAL(wp), DIMENSION(jpi,jpj) :: zwy, zv_trd, zv_frc, zv_spg, zhdiv 162 164 REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhust_e, zhtp2_e 163 165 REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zhvst_e 164 166 REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v ! top/bottom stress at u- & v-points 167 REAL(wp), DIMENSION(jpi,jpj) :: zhU, zhV ! fluxes 165 168 ! 166 169 REAL(wp) :: zwdramp ! local scalar - only used if ln_wd_dl = .True. … … 182 185 zwdramp = r_rn_wdmin1 ! simplest ramp 183 186 ! zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp 184 ! ! reciprocal of baroclinic time step 185 IF( kt == nit000 .AND. neuler == 0 ) THEN ; z2dt_bf = rdt 186 ELSE ; z2dt_bf = 2.0_wp * rdt 187 ENDIF 188 r1_2dt_b = 1.0_wp / z2dt_bf 187 ! ! inverse of baroclinic time step 188 IF( kt == nit000 .AND. neuler == 0 ) THEN ; r1_2dt_b = 1._wp / ( rdt ) 189 ELSE ; r1_2dt_b = 1._wp / ( 2._wp * rdt ) 190 ENDIF 189 191 ! 190 192 ll_init = ln_bt_av ! if no time averaging, then no specific restart … … 210 212 ll_fw_start =.FALSE. 211 213 ENDIF 212 ! 213 ! Set averaging weights and cycle length: 214 ! ! Set averaging weights and cycle length: 214 215 CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) 215 216 ! 216 ENDIF217 !218 IF( ln_isfcav ) THEN ! top+bottom friction (ocean cavities)219 DO jj = 2, jpjm1220 DO ji = fs_2, fs_jpim1 ! vector opt.221 zCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) )222 zCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) )223 END DO224 END DO225 ELSE ! bottom friction only226 DO jj = 2, jpjm1227 DO ji = fs_2, fs_jpim1 ! vector opt.228 zCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )229 zCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )230 END DO231 END DO232 ENDIF233 !234 ! Set arrays to remove/compute coriolis trend.235 ! Do it once at kt=nit000 if volume is fixed, else at each long time step.236 ! Note that these arrays are also used during barotropic loop. These are however frozen237 ! although they should be updated in the variable volume case. Not a big approximation.238 ! To remove this approximation, copy lines below inside barotropic loop239 ! and update depths at T-F points (ht and zhf resp.) at each barotropic time step240 !241 IF( kt == nit000 .OR. .NOT.ln_linssh ) THEN242 !243 SELECT CASE( nvor_scheme )244 CASE( np_EEN ) != EEN scheme using e3f (energy & enstrophy scheme)245 SELECT CASE( nn_een_e3f ) !* ff_f/e3 at F-point246 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)247 DO jj = 1, jpjm1248 DO ji = 1, jpim1249 zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + &250 & ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) * 0.25_wp251 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)252 END DO253 END DO254 CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)255 DO jj = 1, jpjm1256 DO ji = 1, jpim1257 zwz(ji,jj) = ( ht_n (ji ,jj+1) + ht_n (ji+1,jj+1) &258 & + ht_n (ji ,jj ) + ht_n (ji+1,jj ) ) &259 & / ( MAX( 1._wp, ssmask(ji ,jj+1) + ssmask(ji+1,jj+1) &260 & + ssmask(ji ,jj ) + ssmask(ji+1,jj ) ) )261 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj)262 END DO263 END DO264 END SELECT265 CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp )266 !267 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp268 DO jj = 2, jpj269 DO ji = 2, jpi270 ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)271 ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )272 ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)273 ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )274 END DO275 END DO276 !277 CASE( np_EET ) != EEN scheme using e3t (energy conserving scheme)278 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp279 DO jj = 2, jpj280 DO ji = 2, jpi281 z1_ht = ssmask(ji,jj) / ( ht_n(ji,jj) + 1._wp - ssmask(ji,jj) )282 ftne(ji,jj) = ( ff_f(ji-1,jj ) + ff_f(ji ,jj ) + ff_f(ji ,jj-1) ) * z1_ht283 ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) + ff_f(ji ,jj ) ) * z1_ht284 ftse(ji,jj) = ( ff_f(ji ,jj ) + ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht285 ftsw(ji,jj) = ( ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) ) * z1_ht286 END DO287 END DO288 !289 CASE( np_ENE, np_ENS , np_MIX ) != all other schemes (ENE, ENS, MIX) except ENT !290 !291 zwz(:,:) = 0._wp292 zhf(:,:) = 0._wp293 294 !!gm assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed295 !!gm A priori a better value should be something like :296 !!gm zhf(i,j) = masked sum of ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1)297 !!gm divided by the sum of the corresponding mask298 !!gm299 !!300 IF( .NOT.ln_sco ) THEN301 302 !!gm agree the JC comment : this should be done in a much clear way303 304 ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case305 ! Set it to zero for the time being306 ! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level307 ! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth308 ! ENDIF309 ! zhf(:,:) = gdepw_0(:,:,jk+1)310 !311 ELSE312 !313 !zhf(:,:) = hbatf(:,:)314 DO jj = 1, jpjm1315 DO ji = 1, jpim1316 zhf(ji,jj) = ( ht_0 (ji,jj ) + ht_0 (ji+1,jj ) &317 & + ht_0 (ji,jj+1) + ht_0 (ji+1,jj+1) ) &318 & / MAX( ssmask(ji,jj ) + ssmask(ji+1,jj ) &319 & + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp )320 END DO321 END DO322 ENDIF323 !324 DO jj = 1, jpjm1325 zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1))326 END DO327 !328 DO jk = 1, jpkm1329 DO jj = 1, jpjm1330 zhf(:,jj) = zhf(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)331 END DO332 END DO333 CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp )334 ! JC: TBC. hf should be greater than 0335 DO jj = 1, jpj336 DO ji = 1, jpi337 IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj) ! zhf is actually hf here but it saves an array338 END DO339 END DO340 zwz(:,:) = ff_f(:,:) * zwz(:,:)341 END SELECT342 217 ENDIF 343 218 ! … … 348 223 CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) 349 224 ENDIF 225 ! 350 226 351 227 ! ----------------------------------------------------------------------------- … … 354 230 ! 355 231 ! 356 ! !* e3*d/dt(Ua) (Vertically integrated) 357 ! ! -------------------------------------------------- 358 zu_frc(:,:) = 0._wp 359 zv_frc(:,:) = 0._wp 360 ! 361 DO jk = 1, jpkm1 362 zu_frc(:,:) = zu_frc(:,:) + e3u_n(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 363 zv_frc(:,:) = zv_frc(:,:) + e3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 364 END DO 365 ! 366 zu_frc(:,:) = zu_frc(:,:) * r1_hu_n(:,:) 367 zv_frc(:,:) = zv_frc(:,:) * r1_hv_n(:,:) 368 ! 369 ! 370 ! !* baroclinic momentum trend (remove the vertical mean trend) 371 DO jk = 1, jpkm1 ! ----------------------------------------------------------- 372 DO jj = 2, jpjm1 373 DO ji = fs_2, fs_jpim1 ! vector opt. 374 ua(ji,jj,jk) = ua(ji,jj,jk) - zu_frc(ji,jj) * umask(ji,jj,jk) 375 va(ji,jj,jk) = va(ji,jj,jk) - zv_frc(ji,jj) * vmask(ji,jj,jk) 376 END DO 377 END DO 232 ! != zu_frc = 1/H e3*d/dt(Ua) =! (Vertical mean of Ua, the 3D trends) 233 ! ! --------------------------- ! 234 zu_frc(:,:) = SUM( e3u_n(:,:,:) * ua(:,:,:) * umask(:,:,:) , DIM=3 ) * r1_hu_n(:,:) 235 zv_frc(:,:) = SUM( e3v_n(:,:,:) * va(:,:,:) * vmask(:,:,:) , DIM=3 ) * r1_hv_n(:,:) 236 ! 237 ! 238 ! != Ua => baroclinic trend =! (remove its vertical mean) 239 DO jk = 1, jpkm1 ! ------------------------ ! 240 ua(:,:,jk) = ( ua(:,:,jk) - zu_frc(:,:) ) * umask(:,:,jk) 241 va(:,:,jk) = ( va(:,:,jk) - zv_frc(:,:) ) * vmask(:,:,jk) 378 242 END DO 379 243 … … 381 245 !!gm Is it correct to do so ? I think so... 382 246 383 384 ! !* barotropic Coriolis trends (vorticity scheme dependent) 385 ! ! -------------------------------------------------------- 386 ! 387 zwx(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:) ! now fluxes 388 zwy(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:) 389 ! 390 SELECT CASE( nvor_scheme ) 391 CASE( np_ENT ) ! enstrophy conserving scheme (f-point) 392 DO jj = 2, jpjm1 393 DO ji = 2, jpim1 ! vector opt. 394 zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * r1_hu_n(ji,jj) & 395 & * ( e1e2t(ji+1,jj)*ht_n(ji+1,jj)*ff_t(ji+1,jj) * ( vn_b(ji+1,jj) + vn_b(ji+1,jj-1) ) & 396 & + e1e2t(ji ,jj)*ht_n(ji ,jj)*ff_t(ji ,jj) * ( vn_b(ji ,jj) + vn_b(ji ,jj-1) ) ) 397 ! 398 zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * r1_hv_n(ji,jj) & 399 & * ( e1e2t(ji,jj+1)*ht_n(ji,jj+1)*ff_t(ji,jj+1) * ( un_b(ji,jj+1) + un_b(ji-1,jj+1) ) & 400 & + e1e2t(ji,jj )*ht_n(ji,jj )*ff_t(ji,jj ) * ( un_b(ji,jj ) + un_b(ji-1,jj ) ) ) 401 END DO 402 END DO 403 ! 404 CASE( np_ENE , np_MIX ) ! energy conserving scheme (t-point) ENE or MIX 405 DO jj = 2, jpjm1 406 DO ji = fs_2, fs_jpim1 ! vector opt. 407 zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj) 408 zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 409 zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj) 410 zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 411 ! energy conserving formulation for planetary vorticity term 412 zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 413 zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 414 END DO 415 END DO 416 ! 417 CASE( np_ENS ) ! enstrophy conserving scheme (f-point) 418 DO jj = 2, jpjm1 419 DO ji = fs_2, fs_jpim1 ! vector opt. 420 zy1 = r1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 421 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 422 zx1 = - r1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 423 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 424 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 425 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 426 END DO 427 END DO 428 ! 429 CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f) 430 DO jj = 2, jpjm1 431 DO ji = fs_2, fs_jpim1 ! vector opt. 432 zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) & 433 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 434 & + ftse(ji,jj ) * zwy(ji ,jj-1) & 435 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 436 zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 437 & + ftse(ji,jj+1) * zwx(ji ,jj+1) & 438 & + ftnw(ji,jj ) * zwx(ji-1,jj ) & 439 & + ftne(ji,jj ) * zwx(ji ,jj ) ) 440 END DO 441 END DO 442 ! 443 END SELECT 444 ! 445 ! !* Right-Hand-Side of the barotropic momentum equation 446 ! ! ---------------------------------------------------- 447 IF( .NOT.ln_linssh ) THEN ! Variable volume : remove surface pressure gradient 448 IF( ln_wd_il ) THEN ! Calculating and applying W/D gravity filters 247 ! != remove 2D Coriolis and pressure gradient trends =! 248 ! ! ------------------------------------------------- ! 249 ! 250 IF( kt == nit000 .OR. .NOT. ln_linssh ) CALL dyn_cor_2D_init ! Set zwz, the barotropic Coriolis force coefficient 251 ! ! recompute zwz = f/depth at every time step for (.NOT.ln_linssh) as the water colomn height changes 252 ! 253 ! !* 2D Coriolis trends 254 zhU(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:) ! now fluxes 255 zhV(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:) ! NB: FULL domain : put a value in last row and column 256 ! 257 CALL dyn_cor_2d( ht_n, hu_n, hv_n, un_b, vn_b, zhU, zhV, & ! <<== in 258 & zu_trd, zv_trd ) ! ==>> out 259 ! 260 IF( .NOT.ln_linssh ) THEN !* surface pressure gradient (variable volume only) 261 ! 262 IF( ln_wd_il ) THEN ! W/D : limiter applied to spgspg 263 CALL wad_spg( sshn, zcpx, zcpy ) ! Calculating W/D gravity filters, zcpx and zcpy 449 264 DO jj = 2, jpjm1 450 DO ji = 2, jpim1 451 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 452 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 453 & MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji+1,jj) + ht_0(ji+1,jj) ) & 454 & > rn_wdmin1 + rn_wdmin2 455 ll_tmp2 = ( ABS( sshn(ji+1,jj) - sshn(ji ,jj)) > 1.E-12 ).AND.( & 456 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 457 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 458 IF(ll_tmp1) THEN 459 zcpx(ji,jj) = 1.0_wp 460 ELSEIF(ll_tmp2) THEN 461 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 462 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 463 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 464 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 465 ELSE 466 zcpx(ji,jj) = 0._wp 467 ENDIF 468 ! 469 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 470 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 471 & MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji,jj+1) + ht_0(ji,jj+1) ) & 472 & > rn_wdmin1 + rn_wdmin2 473 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1)) > 1.E-12 ).AND.( & 474 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 475 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 476 477 IF(ll_tmp1) THEN 478 zcpy(ji,jj) = 1.0_wp 479 ELSE IF(ll_tmp2) THEN 480 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 481 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 482 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 483 zcpy(ji,jj) = MAX( 0._wp , MIN( zcpy(ji,jj) , 1.0_wp ) ) 484 ELSE 485 zcpy(ji,jj) = 0._wp 486 ENDIF 487 END DO 488 END DO 489 ! 490 DO jj = 2, jpjm1 491 DO ji = 2, jpim1 265 DO ji = 2, jpim1 ! SPG with the application of W/D gravity filters 492 266 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) & 493 267 & * r1_e1u(ji,jj) * zcpx(ji,jj) * wdrampu(ji,jj) !jth … … 496 270 END DO 497 271 END DO 498 ! 499 ELSE 500 ! 272 ELSE ! now suface pressure gradient 501 273 DO jj = 2, jpjm1 502 274 DO ji = fs_2, fs_jpim1 ! vector opt. … … 516 288 END DO 517 289 ! 518 ! ! Add bottom stress contribution from baroclinic velocities: 519 IF (ln_bt_fw) THEN 520 DO jj = 2, jpjm1 521 DO ji = fs_2, fs_jpim1 ! vector opt. 522 ikbu = mbku(ji,jj) 523 ikbv = mbkv(ji,jj) 524 zwx(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) ! NOW bottom baroclinic velocities 525 zwy(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj) 526 END DO 527 END DO 528 ELSE 529 DO jj = 2, jpjm1 530 DO ji = fs_2, fs_jpim1 ! vector opt. 531 ikbu = mbku(ji,jj) 532 ikbv = mbkv(ji,jj) 533 zwx(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) ! BEFORE bottom baroclinic velocities 534 zwy(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj) 535 END DO 536 END DO 537 ENDIF 538 ! 539 ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 540 IF( ln_wd_il ) THEN 541 zztmp = -1._wp / rdtbt 542 DO jj = 2, jpjm1 543 DO ji = fs_2, fs_jpim1 ! vector opt. 544 zu_frc(ji,jj) = zu_frc(ji,jj) + & 545 & MAX(r1_hu_n(ji,jj) * r1_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ), zztmp ) * zwx(ji,jj) * wdrampu(ji,jj) 546 zv_frc(ji,jj) = zv_frc(ji,jj) + & 547 & MAX(r1_hv_n(ji,jj) * r1_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ), zztmp ) * zwy(ji,jj) * wdrampv(ji,jj) 548 END DO 549 END DO 550 ELSE 551 DO jj = 2, jpjm1 552 DO ji = fs_2, fs_jpim1 ! vector opt. 553 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * r1_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zwx(ji,jj) 554 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * r1_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zwy(ji,jj) 555 END DO 556 END DO 557 END IF 558 ! 559 IF( ln_isfcav ) THEN ! Add TOP stress contribution from baroclinic velocities: 560 IF( ln_bt_fw ) THEN 561 DO jj = 2, jpjm1 290 ! != Add bottom stress contribution from baroclinic velocities =! 291 ! ! ----------------------------------------------------------- ! 292 CALL drg_init( zu_frc, zv_frc, zCdU_u, zCdU_v ) ! also provide the barotropic drag coefficients 293 ! 294 ! != Add atmospheric pressure forcing =! 295 ! ! ---------------------------------- ! 296 IF( ln_apr_dyn ) THEN 297 IF( ln_bt_fw ) THEN ! FORWARD integration: use kt+1/2 pressure (NOW+1/2) 298 DO jj = 2, jpjm1 562 299 DO ji = fs_2, fs_jpim1 ! vector opt. 563 iktu = miku(ji,jj) 564 iktv = mikv(ji,jj) 565 zwx(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) ! NOW top baroclinic velocities 566 zwy(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj) 567 END DO 568 END DO 569 ELSE 570 DO jj = 2, jpjm1 300 zu_frc(ji,jj) = zu_frc(ji,jj) + grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj) 301 zv_frc(ji,jj) = zv_frc(ji,jj) + grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj) 302 END DO 303 END DO 304 ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 pressure (NOW) 305 zztmp = grav * r1_2 306 DO jj = 2, jpjm1 571 307 DO ji = fs_2, fs_jpim1 ! vector opt. 572 iktu = miku(ji,jj) 573 iktv = mikv(ji,jj) 574 zwx(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) ! BEFORE top baroclinic velocities 575 zwy(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj) 576 END DO 577 END DO 578 ENDIF 579 ! 580 ! Note that the "unclipped" top friction parameter is used even with explicit drag 581 DO jj = 2, jpjm1 582 DO ji = fs_2, fs_jpim1 ! vector opt. 583 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * r1_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zwx(ji,jj) 584 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * r1_2 * ( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zwy(ji,jj) 585 END DO 586 END DO 587 ENDIF 588 ! 308 zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) & 309 & + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj) 310 zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) & 311 & + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj) 312 END DO 313 END DO 314 ENDIF 315 ENDIF 316 ! 317 ! != Add atmospheric pressure forcing =! 318 ! ! ---------------------------------- ! 589 319 IF( ln_bt_fw ) THEN ! Add wind forcing 590 320 DO jj = 2, jpjm1 … … 604 334 ENDIF 605 335 ! 606 IF( ln_apr_dyn ) THEN ! Add atm pressure forcing 607 IF( ln_bt_fw ) THEN 608 DO jj = 2, jpjm1 609 DO ji = fs_2, fs_jpim1 ! vector opt. 610 zu_spg = grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj) 611 zv_spg = grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj) 612 zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg 613 zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg 614 END DO 615 END DO 616 ELSE 617 zztmp = grav * r1_2 618 DO jj = 2, jpjm1 619 DO ji = fs_2, fs_jpim1 ! vector opt. 620 zu_spg = zztmp * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) & 621 & + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj) 622 zv_spg = zztmp * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) & 623 & + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj) 624 zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg 625 zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg 626 END DO 627 END DO 628 ENDIF 629 ENDIF 630 ! !* Right-Hand-Side of the barotropic ssh equation 631 ! ! ----------------------------------------------- 632 ! ! Surface net water flux and rivers 633 IF (ln_bt_fw) THEN 634 zssh_frc(:,:) = r1_rau0 * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) 635 ELSE 336 ! !----------------! 337 ! !== sssh_frc ==! Right-Hand-Side of the barotropic ssh equation (over the FULL domain) 338 ! !----------------! 339 ! != Net water flux forcing applied to a water column =! 340 ! ! --------------------------------------------------- ! 341 IF (ln_bt_fw) THEN ! FORWARD integration: use kt+1/2 fluxes (NOW+1/2) 342 zssh_frc(:,:) = r1_rau0 * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) 343 ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 fluxes (NOW) 636 344 zztmp = r1_rau0 * r1_2 637 zssh_frc(:,:) = zztmp * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) & 638 & + fwfisf(:,:) + fwfisf_b(:,:) ) 639 ENDIF 640 ! 641 IF( ln_sdw ) THEN ! Stokes drift divergence added if necessary 345 zssh_frc(:,:) = zztmp * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) + fwfisf(:,:) + fwfisf_b(:,:) ) 346 ENDIF 347 ! != Add Stokes drift divergence =! (if exist) 348 IF( ln_sdw ) THEN ! ----------------------------- ! 642 349 zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:) 643 350 ENDIF 644 351 ! 645 352 #if defined key_asminc 646 ! ! Include the IAU weighted SSH increment 353 ! != Add the IAU weighted SSH increment =! 354 ! ! ------------------------------------ ! 647 355 IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN 648 356 zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:) 649 357 ENDIF 650 358 #endif 651 ! ! *Fill boundary data arrays for AGRIF359 ! != Fill boundary data arrays for AGRIF 652 360 ! ! ------------------------------------ 653 361 #if defined key_agrif … … 671 379 vb_e (:,:) = 0._wp 672 380 ENDIF 673 381 ! 382 IF( ln_linssh ) THEN ! mid-step ocean depth is fixed (hup2_e=hu_n=hu_0) 383 zhup2_e(:,:) = hu_n(:,:) 384 zhvp2_e(:,:) = hv_n(:,:) 385 zhtp2_e(:,:) = ht_n(:,:) 386 ENDIF 674 387 ! 675 388 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields … … 693 406 ENDIF 694 407 ! 695 !696 !697 408 ! Initialize sums: 698 409 ua_b (:,:) = 0._wp ! After barotropic velocities (or transport if flux form) … … 714 425 ! 715 426 l_full_nf_update = jn == icycle ! false: disable full North fold update (performances) for jn = 1 to icycle-1 716 ! ! ------------------ 717 ! !* Update the forcing (BDY and tides) 718 ! ! ------------------ 719 ! Update only tidal forcing at open boundaries 427 ! 428 ! !== Update the forcing ==! (BDY and tides) 429 ! 720 430 IF( ln_bdy .AND. ln_tide ) CALL bdy_dta_tides( kt, kit=jn, kt_offset= noffset+1 ) 721 431 IF( ln_tide_pot .AND. ln_tide ) CALL upd_tide ( kt, kit=jn, kt_offset= noffset ) 722 432 ! 723 ! Set extrapolation coefficients for predictor step: 433 ! !== extrapolation at mid-step ==! (jn+1/2) 434 ! 435 ! !* Set extrapolation coefficients for predictor step: 724 436 IF ((jn<3).AND.ll_init) THEN ! Forward 725 437 za1 = 1._wp … … 731 443 za3 = 0.281105_wp ! za3 = bet 732 444 ENDIF 733 734 ! Extrapolate barotropic velocities at step jit+0.5: 445 ! 446 ! !* Extrapolate barotropic velocities at mid-step (jn+1/2) 447 !-- m+1/2 m m-1 m-2 --! 448 !-- u = (3/2+beta) u -(1/2+2beta) u + beta u --! 449 !-------------------------------------------------------------------------! 735 450 ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:) 736 451 va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:) … … 739 454 ! ! ------------------ 740 455 ! Extrapolate Sea Level at step jit+0.5: 456 !-- m+1/2 m m-1 m-2 --! 457 !-- ssh = (3/2+beta) ssh -(1/2+2beta) ssh + beta ssh --! 458 !--------------------------------------------------------------------------------! 741 459 zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 742 460 743 ! set wetting & drying mask at tracer points for this barotropic sub-step 744 IF ( ln_wd_dl ) THEN 745 ! 746 IF ( ln_wd_dl_rmp ) THEN 747 DO jj = 1, jpj 748 DO ji = 1, jpi ! vector opt. 749 IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN 750 ! IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN 751 ztwdmask(ji,jj) = 1._wp 752 ELSE IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN 753 ztwdmask(ji,jj) = (tanh(50._wp*( ( zsshp2_e(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )*r_rn_wdmin1)) ) 754 ELSE 755 ztwdmask(ji,jj) = 0._wp 756 END IF 757 END DO 758 END DO 759 ELSE 760 DO jj = 1, jpj 761 DO ji = 1, jpi ! vector opt. 762 IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN 763 ztwdmask(ji,jj) = 1._wp 764 ELSE 765 ztwdmask(ji,jj) = 0._wp 766 ENDIF 767 END DO 768 END DO 769 ENDIF 770 ! 771 ENDIF 461 ! set wetting & drying mask at tracer points for this barotropic mid-step 462 IF( ln_wd_dl ) CALL wad_tmsk( zsshp2_e, ztwdmask ) 772 463 ! 773 464 DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points … … 786 477 zhvp2_e(:,:) = hv_0(:,:) + zwy(:,:) 787 478 zhtp2_e(:,:) = ht_0(:,:) + zsshp2_e(:,:) 788 ELSE 789 zhup2_e(:,:) = hu_n(:,:) 790 zhvp2_e(:,:) = hv_n(:,:) 791 zhtp2_e(:,:) = ht_n(:,:) 792 ENDIF 793 ! !* after ssh 794 ! ! ----------- 795 ! 796 ! Enforce volume conservation at open boundaries: 479 ! 480 ENDIF 481 ! 482 ! !== after SSH ==! (jn+1) 483 ! 484 ! ! update (ua_e,va_e) to enforce volume conservation at open boundaries 485 ! ! values of zhup2_e and zhvp2_e on the halo are not needed in bdy_vol2d 797 486 IF( ln_bdy .AND. ln_vol ) CALL bdy_vol2d( kt, jn, ua_e, va_e, zhup2_e, zhvp2_e ) 798 487 ! 799 z wx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5800 z wy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:)488 zhU(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5 489 zhV(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 801 490 ! 802 491 #if defined key_agrif … … 805 494 IF((nbondi == -1).OR.(nbondi == 2)) THEN 806 495 DO jj = 1, jpj 807 z wx(2:nbghostcells+1,jj) = ubdy_w(1:nbghostcells,jj) * e2u(2:nbghostcells+1,jj)808 z wy(2:nbghostcells+1,jj) = vbdy_w(1:nbghostcells,jj) * e1v(2:nbghostcells+1,jj)496 zhU(2:nbghostcells+1,jj) = ubdy_w(1:nbghostcells,jj) * e2u(2:nbghostcells+1,jj) 497 zhV(2:nbghostcells+1,jj) = vbdy_w(1:nbghostcells,jj) * e1v(2:nbghostcells+1,jj) 809 498 END DO 810 499 ENDIF 811 500 IF((nbondi == 1).OR.(nbondi == 2)) THEN 812 501 DO jj=1,jpj 813 z wx(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(1:nbghostcells,jj) * e2u(nlci-nbghostcells-1:nlci-2,jj)814 z wy(nlci-nbghostcells :nlci-1,jj) = vbdy_e(1:nbghostcells,jj) * e1v(nlci-nbghostcells :nlci-1,jj)502 zhU(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(1:nbghostcells,jj) * e2u(nlci-nbghostcells-1:nlci-2,jj) 503 zhV(nlci-nbghostcells :nlci-1,jj) = vbdy_e(1:nbghostcells,jj) * e1v(nlci-nbghostcells :nlci-1,jj) 815 504 END DO 816 505 ENDIF 817 506 IF((nbondj == -1).OR.(nbondj == 2)) THEN 818 507 DO ji=1,jpi 819 z wy(ji,2:nbghostcells+1) = vbdy_s(ji,1:nbghostcells) * e1v(ji,2:nbghostcells+1)820 z wx(ji,2:nbghostcells+1) = ubdy_s(ji,1:nbghostcells) * e2u(ji,2:nbghostcells+1)508 zhV(ji,2:nbghostcells+1) = vbdy_s(ji,1:nbghostcells) * e1v(ji,2:nbghostcells+1) 509 zhU(ji,2:nbghostcells+1) = ubdy_s(ji,1:nbghostcells) * e2u(ji,2:nbghostcells+1) 821 510 END DO 822 511 ENDIF 823 512 IF((nbondj == 1).OR.(nbondj == 2)) THEN 824 513 DO ji=1,jpi 825 z wy(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji,1:nbghostcells) * e1v(ji,nlcj-nbghostcells-1:nlcj-2)826 z wx(ji,nlcj-nbghostcells :nlcj-1) = ubdy_n(ji,1:nbghostcells) * e2u(ji,nlcj-nbghostcells :nlcj-1)514 zhV(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji,1:nbghostcells) * e1v(ji,nlcj-nbghostcells-1:nlcj-2) 515 zhU(ji,nlcj-nbghostcells :nlcj-1) = ubdy_n(ji,1:nbghostcells) * e2u(ji,nlcj-nbghostcells :nlcj-1) 827 516 END DO 828 517 ENDIF 829 518 ENDIF 830 519 #endif 831 IF( ln_wd_il ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 832 833 IF ( ln_wd_dl ) THEN 834 ! 835 ! un_e and vn_e are set to zero at faces where the direction of the flow is from dry cells 836 ! 837 DO jj = 1, jpjm1 838 DO ji = 1, jpim1 839 IF ( zwx(ji,jj) > 0.0 ) THEN 840 zuwdmask(ji, jj) = ztwdmask(ji ,jj) 841 ELSE 842 zuwdmask(ji, jj) = ztwdmask(ji+1,jj) 843 END IF 844 zwx(ji, jj) = zuwdmask(ji,jj)*zwx(ji, jj) 845 un_e(ji,jj) = zuwdmask(ji,jj)*un_e(ji,jj) 846 847 IF ( zwy(ji,jj) > 0.0 ) THEN 848 zvwdmask(ji, jj) = ztwdmask(ji, jj ) 849 ELSE 850 zvwdmask(ji, jj) = ztwdmask(ji, jj+1) 851 END IF 852 zwy(ji, jj) = zvwdmask(ji,jj)*zwy(ji,jj) 853 vn_e(ji,jj) = zvwdmask(ji,jj)*vn_e(ji,jj) 854 END DO 855 END DO 520 IF( ln_wd_il ) CALL wad_lmt_bt(zhU, zhV, sshn_e, zssh_frc, rdtbt) 521 522 IF( ln_wd_dl ) THEN ! un_e and vn_e are set to zero at faces where 523 ! ! the direction of the flow is from dry cells 524 CALL wad_Umsk( ztwdmask, zhU, zhV, un_e, vn_e, zuwdmask, zvwdmask ) 856 525 ! 857 526 ENDIF 858 859 ! Sum over sub-time-steps to compute advective velocities 860 za2 = wgtbtp2(jn) 861 un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:) 862 vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:) 863 864 ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc = True) 527 ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True) 865 528 IF ( ln_wd_dl_bc ) THEN 866 529 zuwdav2(:,:) = zuwdav2(:,:) + za2 * zuwdmask(:,:) 867 530 zvwdav2(:,:) = zvwdav2(:,:) + za2 * zvwdmask(:,:) 868 531 END IF 532 533 ! Sum over sub-time-steps to compute advective velocities 534 ! 535 za2 = wgtbtp2(jn) 536 un_adv(:,:) = un_adv(:,:) + za2 * zhU(:,:) * r1_e2u(:,:) 537 vn_adv(:,:) = vn_adv(:,:) + za2 * zhV(:,:) * r1_e1v(:,:) 869 538 870 539 ! Set next sea level: 871 540 DO jj = 2, jpjm1 872 541 DO ji = fs_2, fs_jpim1 ! vector opt. 873 zhdiv(ji,jj) = ( zwx(ji,jj) - zwx(ji-1,jj) & 874 & + zwy(ji,jj) - zwy(ji,jj-1) ) * r1_e1e2t(ji,jj) 875 END DO 876 END DO 542 zhdiv(ji,jj) = ( zhU(ji,jj) - zhU(ji-1,jj) & 543 & + zhV(ji,jj) - zhV(ji,jj-1) ) * r1_e1e2t(ji,jj) 544 END DO 545 END DO 546 ! Compute Sea Level at step jit+1 547 !-- m+1 m m+1/2 --! 548 !-- ssh = ssh - delta_t' * [ frc + div( flux ) ] --! 549 !-------------------------------------------------------------------------! 877 550 ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * ssmask(:,:) 878 551 … … 899 572 CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) 900 573 ENDIF 901 ! 902 ! Half-step back interpolation of SSH for surface pressure computation: 903 !---------------------------------------------------------------------- 904 IF ((jn==1).AND.ll_init) THEN 905 za0=1._wp ! Forward-backward 906 za1=0._wp 907 za2=0._wp 908 za3=0._wp 909 ELSEIF ((jn==2).AND.ll_init) THEN ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12 910 za0= 1.0833333333333_wp ! za0 = 1-gam-eps 911 za1=-0.1666666666666_wp ! za1 = gam 912 za2= 0.0833333333333_wp ! za2 = eps 913 za3= 0._wp 914 ELSE ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880 915 IF (rn_bt_alpha==0._wp) THEN 916 za0=0.614_wp ! za0 = 1/2 + gam + 2*eps 917 za1=0.285_wp ! za1 = 1/2 - 2*gam - 3*eps 918 za2=0.088_wp ! za2 = gam 919 za3=0.013_wp ! za3 = eps 920 ELSE 921 zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha 922 zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha 923 za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon 924 za1 = 1._wp - za0 - zgamma - zepsilon 925 za2 = zgamma 926 za3 = zepsilon 927 ENDIF 928 ENDIF 929 ! 574 ! 575 ! Half-step back interpolation of SSH for surface pressure computation at step jit+1/2 576 !-- m+1/2 m+1 m m-1 m-2 --! 577 !-- ssh' = za0 * ssh + za1 * ssh + za2 * ssh + za3 * ssh --! 578 !------------------------------------------------------------------------------------------! 579 CALL ts_bck_interp( jn, ll_init, za0, za1, za2, za3 ) ! coeficients of the interpolation 930 580 zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) & 931 581 & + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 932 933 IF( ln_wd_il ) THEN ! Calculating and applying W/D gravity filters 934 DO jj = 2, jpjm1 935 DO ji = 2, jpim1 936 ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > & 937 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 938 & MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) ) & 939 & > rn_wdmin1 + rn_wdmin2 940 ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji+1,jj)) > 1.E-12 ).AND.( & 941 & MAX( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > & 942 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 943 944 IF(ll_tmp1) THEN 945 zcpx(ji,jj) = 1.0_wp 946 ELSE IF(ll_tmp2) THEN 947 ! no worries about zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj) = 0, it won't happen ! here 948 zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) - zsshp2_e(ji,jj) - ht_0(ji,jj)) & 949 & / (zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj)) ) 950 ELSE 951 zcpx(ji,jj) = 0._wp 952 ENDIF 953 ! 954 ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > & 955 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 956 & MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) ) & 957 & > rn_wdmin1 + rn_wdmin2 958 ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji,jj+1)) > 1.E-12 ).AND.( & 959 & MAX( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > & 960 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 961 962 IF(ll_tmp1) THEN 963 zcpy(ji,jj) = 1.0_wp 964 ELSEIF(ll_tmp2) THEN 965 ! no worries about zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj ) = 0, it won't happen ! here 966 zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) - zsshp2_e(ji,jj) - ht_0(ji,jj)) & 967 & / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj )) ) 968 ELSE 969 zcpy(ji,jj) = 0._wp 970 ENDIF 971 END DO 972 END DO 973 ENDIF 974 ! 975 ! Compute associated depths at U and V points: 976 IF( .NOT.ln_linssh .AND. .NOT.ln_dynadv_vec ) THEN !* Vector form 977 ! 978 DO jj = 2, jpjm1 979 DO ji = 2, jpim1 980 zx1 = r1_2 * ssumask(ji ,jj) * r1_e1e2u(ji ,jj) & 981 & * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj) & 982 & + e1e2t(ji+1,jj ) * zsshp2_e(ji+1,jj ) ) 983 zy1 = r1_2 * ssvmask(ji ,jj) * r1_e1e2v(ji ,jj ) & 984 & * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj ) & 985 & + e1e2t(ji ,jj+1) * zsshp2_e(ji ,jj+1) ) 986 zhust_e(ji,jj) = hu_0(ji,jj) + zx1 987 zhvst_e(ji,jj) = hv_0(ji,jj) + zy1 988 END DO 989 END DO 990 ! 582 ! 583 ! ! Surface pressure gradient 584 zldg = ( 1._wp - rn_scal_load ) * grav ! local factor 585 DO jj = 2, jpjm1 586 DO ji = 2, jpim1 587 zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 588 zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 589 END DO 590 END DO 591 IF( ln_wd_il ) THEN ! W/D : gravity filters applied on pressure gradient 592 CALL wad_spg( zsshp2_e, zcpx, zcpy ) ! Calculating W/D gravity filters 593 zu_spg(2:jpim1,2:jpjm1) = zu_spg(2:jpim1,2:jpjm1) * zcpx(2:jpim1,2:jpjm1) 594 zv_spg(2:jpim1,2:jpjm1) = zv_spg(2:jpim1,2:jpjm1) * zcpy(2:jpim1,2:jpjm1) 991 595 ENDIF 992 596 ! … … 994 598 ! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated 995 599 ! at each time step. We however keep them constant here for optimization. 996 ! Recall that zwx and zwy arrays hold fluxes at this stage: 997 ! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5 998 ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 999 ! 1000 SELECT CASE( nvor_scheme ) 1001 CASE( np_ENT ) ! energy conserving scheme (t-point) 1002 DO jj = 2, jpjm1 1003 DO ji = 2, jpim1 ! vector opt. 1004 1005 z1_hu = ssumask(ji,jj) / ( zhup2_e(ji,jj) + 1._wp - ssumask(ji,jj) ) 1006 z1_hv = ssvmask(ji,jj) / ( zhvp2_e(ji,jj) + 1._wp - ssvmask(ji,jj) ) 1007 1008 zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu & 1009 & * ( e1e2t(ji+1,jj)*zhtp2_e(ji+1,jj)*ff_t(ji+1,jj) * ( va_e(ji+1,jj) + va_e(ji+1,jj-1) ) & 1010 & + e1e2t(ji ,jj)*zhtp2_e(ji ,jj)*ff_t(ji ,jj) * ( va_e(ji ,jj) + va_e(ji ,jj-1) ) ) 1011 ! 1012 zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv & 1013 & * ( e1e2t(ji,jj+1)*zhtp2_e(ji,jj+1)*ff_t(ji,jj+1) * ( ua_e(ji,jj+1) + ua_e(ji-1,jj+1) ) & 1014 & + e1e2t(ji,jj )*zhtp2_e(ji,jj )*ff_t(ji,jj ) * ( ua_e(ji,jj ) + ua_e(ji-1,jj ) ) ) 1015 END DO 1016 END DO 1017 ! 1018 CASE( np_ENE, np_MIX ) ! energy conserving scheme (f-point) 1019 DO jj = 2, jpjm1 1020 DO ji = fs_2, fs_jpim1 ! vector opt. 1021 zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj) 1022 zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 1023 zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj) 1024 zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 1025 zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 1026 zv_trd(ji,jj) =-r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 1027 END DO 1028 END DO 1029 ! 1030 CASE( np_ENS ) ! enstrophy conserving scheme (f-point) 1031 DO jj = 2, jpjm1 1032 DO ji = fs_2, fs_jpim1 ! vector opt. 1033 zy1 = r1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 1034 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 1035 zx1 = - r1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 1036 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 1037 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 1038 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 1039 END DO 1040 END DO 1041 ! 1042 CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f) 1043 DO jj = 2, jpjm1 1044 DO ji = fs_2, fs_jpim1 ! vector opt. 1045 zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) & 1046 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 1047 & + ftse(ji,jj ) * zwy(ji ,jj-1) & 1048 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 1049 zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 1050 & + ftse(ji,jj+1) * zwx(ji ,jj+1) & 1051 & + ftnw(ji,jj ) * zwx(ji-1,jj ) & 1052 & + ftne(ji,jj ) * zwx(ji ,jj ) ) 1053 END DO 1054 END DO 1055 ! 1056 END SELECT 600 ! Recall that zhU and zhV hold fluxes at jn+0.5 (extrapolated not backward interpolated) 601 CALL dyn_cor_2d( zhtp2_e, zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV, zu_trd, zv_trd ) 1057 602 ! 1058 603 ! Add tidal astronomical forcing if defined … … 1060 605 DO jj = 2, jpjm1 1061 606 DO ji = fs_2, fs_jpim1 ! vector opt. 1062 zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 1063 zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) 1064 zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg 1065 zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg 607 zu_trd(ji,jj) = zu_trd(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 608 zv_trd(ji,jj) = zv_trd(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) 1066 609 END DO 1067 610 END DO … … 1077 620 END DO 1078 621 END DO 1079 ENDIF 1080 ! 1081 ! Surface pressure trend: 1082 IF( ln_wd_il ) THEN 1083 DO jj = 2, jpjm1 1084 DO ji = 2, jpim1 1085 ! Add surface pressure gradient 1086 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 1087 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 1088 zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg * zcpx(ji,jj) 1089 zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg * zcpy(ji,jj) 1090 END DO 1091 END DO 1092 ELSE 1093 DO jj = 2, jpjm1 1094 DO ji = fs_2, fs_jpim1 ! vector opt. 1095 ! Add surface pressure gradient 1096 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 1097 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 1098 zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg 1099 zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg 1100 END DO 1101 END DO 1102 END IF 1103 622 ENDIF 1104 623 ! 1105 624 ! Set next velocities: 625 ! Compute barotropic speeds at step jit+1 (h : total height of the water colomn) 626 !-- VECTOR FORM 627 !-- m+1 m / m+1/2 \ --! 628 !-- u = u + delta_t' * \ grad_x( ssh') - f * k vect u + frc / --! 629 !-- --! 630 !-- FLUX FORM --! 631 !-- m+1 ___1___ / m m / m+1/2 m+1/2 n \ \ --! 632 !-- u = m+1 | h * u + delta_t' * \ grad_x( ssh') - h * f * k vect u + h * frc / | --! 633 !-- h \ / --! 634 !------------------------------------------------------------------------------------------------------------! 1106 635 IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form 1107 636 DO jj = 2, jpjm1 1108 637 DO ji = fs_2, fs_jpim1 ! vector opt. 1109 638 ua_e(ji,jj) = ( un_e(ji,jj) & 1110 & + rdtbt * ( zwx(ji,jj) &639 & + rdtbt * ( zu_spg(ji,jj) & 1111 640 & + zu_trd(ji,jj) & 1112 641 & + zu_frc(ji,jj) ) & … … 1114 643 1115 644 va_e(ji,jj) = ( vn_e(ji,jj) & 1116 & + rdtbt * ( zwy(ji,jj) &645 & + rdtbt * ( zv_spg(ji,jj) & 1117 646 & + zv_trd(ji,jj) & 1118 647 & + zv_frc(ji,jj) ) & 1119 648 & ) * ssvmask(ji,jj) 1120 1121 649 END DO 1122 650 END DO … … 1125 653 DO jj = 2, jpjm1 1126 654 DO ji = fs_2, fs_jpim1 ! vector opt. 1127 1128 zhu ra = hu_0(ji,jj) + zsshu_a(ji,jj)1129 zhvra = hv_0(ji,jj) + zsshv_a(ji,jj)1130 1131 zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj))1132 zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj))1133 1134 ua_e(ji,jj) = ( hu_e(ji,jj) * un_e(ji,jj) &1135 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) &1136 & + zhup2_e(ji,jj) * zu_trd(ji,jj) &1137 & + hu_n(ji,jj) * zu_frc(ji,jj) ) &1138 & ) * zhura1139 1140 va_e(ji,jj) = ( hv_e(ji,jj) * vn_e(ji,jj) &1141 & + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj)&1142 & + zhvp2_e(ji,jj) * zv_trd(ji,jj) &1143 & + hv_n(ji,jj) * zv_frc(ji,jj) ) &1144 & ) * zhvra655 ! ! backward extrapolated depth used in spg terms at jn+1/2 656 zhu_bck = hu_0(ji,jj) + r1_2*r1_e1e2u(ji,jj) * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 657 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj) 658 zhv_bck = hv_0(ji,jj) + r1_2*r1_e1e2v(ji,jj) * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 659 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj) 660 ! ! inverse depth at jn+1 661 z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) ) 662 z1_hv = ssvmask(ji,jj) / ( hu_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj) ) 663 ! 664 ua_e(ji,jj) = ( hu_e (ji,jj) * un_e (ji,jj) & 665 & + rdtbt * ( zhu_bck * zu_spg (ji,jj) & ! 666 & + zhup2_e(ji,jj) * zu_trd (ji,jj) & ! 667 & + hu_n (ji,jj) * zu_frc (ji,jj) ) ) * z1_hu 668 ! 669 va_e(ji,jj) = ( hv_e (ji,jj) * vn_e (ji,jj) & 670 & + rdtbt * ( zhv_bck * zv_spg (ji,jj) & ! 671 & + zhvp2_e(ji,jj) * zv_trd (ji,jj) & ! 672 & + hv_n (ji,jj) * zv_frc (ji,jj) ) ) * z1_hu 1145 673 END DO 1146 674 END DO … … 1213 741 ! Set advection velocity correction: 1214 742 IF (ln_bt_fw) THEN 1215 zwx(:,:) = un_adv(:,:)1216 zwy(:,:) = vn_adv(:,:)1217 743 IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) THEN 1218 un_adv(:,:) = r1_2 * ( ub2_b(:,:) + zwx(:,:) - atfp * un_bf(:,:) ) 1219 vn_adv(:,:) = r1_2 * ( vb2_b(:,:) + zwy(:,:) - atfp * vn_bf(:,:) ) 1220 ! 1221 ! Update corrective fluxes for next time step: 1222 un_bf(:,:) = atfp * un_bf(:,:) + (zwx(:,:) - ub2_b(:,:)) 1223 vn_bf(:,:) = atfp * vn_bf(:,:) + (zwy(:,:) - vb2_b(:,:)) 744 DO jj = 1, jpj 745 DO ji = 1, jpi 746 zun_save = un_adv(ji,jj) 747 zvn_save = vn_adv(ji,jj) 748 ! ! apply the previously computed correction 749 un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - atfp * un_bf(ji,jj) ) 750 vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - atfp * vn_bf(ji,jj) ) 751 ! ! Update corrective fluxes for next time step 752 un_bf(ji,jj) = atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) ) 753 vn_bf(ji,jj) = atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) ) 754 ! ! Save integrated transport for next computation 755 ub2_b(ji,jj) = zun_save 756 vb2_b(ji,jj) = zvn_save 757 END DO 758 END DO 1224 759 ELSE 1225 un_bf(:,:) = 0._wp 1226 vn_bf(:,:) = 0._wp 1227 END IF 1228 ! Save integrated transport for next computation 1229 ub2_b(:,:) = zwx(:,:) 1230 vb2_b(:,:) = zwy(:,:) 760 un_bf(:,:) = 0._wp ! corrective fluxes for next time step set to zero 761 vn_bf(:,:) = 0._wp 762 ub2_b(:,:) = un_adv(:,:) ! Save integrated transport for next computation 763 vb2_b(:,:) = vn_adv(:,:) 764 END IF 1231 765 ENDIF 1232 766 … … 1473 1007 REAL(wp) :: zxr2, zyr2, zcmax ! local scalar 1474 1008 REAL(wp), DIMENSION(jpi,jpj) :: zcu 1009 INTEGER :: inum 1475 1010 !!---------------------------------------------------------------------- 1476 1011 ! … … 1579 1114 END SUBROUTINE dyn_spg_ts_init 1580 1115 1116 1117 SUBROUTINE dyn_cor_2d_init 1118 !!--------------------------------------------------------------------- 1119 !! *** ROUTINE dyn_cor_2d_init *** 1120 !! 1121 !! ** Purpose : Set time splitting options 1122 !! Set arrays to remove/compute coriolis trend. 1123 !! Do it once during initialization if volume is fixed, else at each long time step. 1124 !! Note that these arrays are also used during barotropic loop. These are however frozen 1125 !! although they should be updated in the variable volume case. Not a big approximation. 1126 !! To remove this approximation, copy lines below inside barotropic loop 1127 !! and update depths at T-F points (ht and zhf resp.) at each barotropic time step 1128 !! 1129 !! Compute zwz = f / ( height of the water colomn ) 1130 !!---------------------------------------------------------------------- 1131 INTEGER :: ji ,jj, jk ! dummy loop indices 1132 REAL(wp) :: z1_ht 1133 REAL(wp), DIMENSION(jpi,jpj) :: zhf 1134 !!---------------------------------------------------------------------- 1135 ! 1136 ! 1137 SELECT CASE( nvor_scheme ) 1138 CASE( np_EEN ) != EEN scheme using e3f (energy & enstrophy scheme) 1139 SELECT CASE( nn_een_e3f ) !* ff_f/e3 at F-point 1140 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 1141 DO jj = 1, jpjm1 1142 DO ji = 1, jpim1 1143 zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + & 1144 & ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) * 0.25_wp 1145 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 1146 END DO 1147 END DO 1148 CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask) 1149 DO jj = 1, jpjm1 1150 DO ji = 1, jpim1 1151 zwz(ji,jj) = ( ht_n (ji ,jj+1) + ht_n (ji+1,jj+1) & 1152 & + ht_n (ji ,jj ) + ht_n (ji+1,jj ) ) & 1153 & / ( MAX( 1._wp, ssmask(ji ,jj+1) + ssmask(ji+1,jj+1) & 1154 & + ssmask(ji ,jj ) + ssmask(ji+1,jj ) ) ) 1155 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 1156 END DO 1157 END DO 1158 END SELECT 1159 CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp ) 1160 ! 1161 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 1162 DO jj = 2, jpj 1163 DO ji = 2, jpi 1164 ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 1165 ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) 1166 ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) 1167 ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) 1168 END DO 1169 END DO 1170 ! 1171 CASE( np_EET ) != EEN scheme using e3t (energy conserving scheme) 1172 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 1173 DO jj = 2, jpj 1174 DO ji = 2, jpi 1175 z1_ht = ssmask(ji,jj) / ( ht_n(ji,jj) + 1._wp - ssmask(ji,jj) ) 1176 ftne(ji,jj) = ( ff_f(ji-1,jj ) + ff_f(ji ,jj ) + ff_f(ji ,jj-1) ) * z1_ht 1177 ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) + ff_f(ji ,jj ) ) * z1_ht 1178 ftse(ji,jj) = ( ff_f(ji ,jj ) + ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht 1179 ftsw(ji,jj) = ( ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) ) * z1_ht 1180 END DO 1181 END DO 1182 ! 1183 CASE( np_ENE, np_ENS , np_MIX ) != all other schemes (ENE, ENS, MIX) except ENT ! 1184 ! 1185 zwz(:,:) = 0._wp 1186 zhf(:,:) = 0._wp 1187 1188 !!gm assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed 1189 !!gm A priori a better value should be something like : 1190 !!gm zhf(i,j) = masked sum of ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1) 1191 !!gm divided by the sum of the corresponding mask 1192 !!gm 1193 !! 1194 IF( .NOT.ln_sco ) THEN 1195 1196 !!gm agree the JC comment : this should be done in a much clear way 1197 1198 ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case 1199 ! Set it to zero for the time being 1200 ! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level 1201 ! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth 1202 ! ENDIF 1203 ! zhf(:,:) = gdepw_0(:,:,jk+1) 1204 ! 1205 ELSE 1206 ! 1207 !zhf(:,:) = hbatf(:,:) 1208 DO jj = 1, jpjm1 1209 DO ji = 1, jpim1 1210 zhf(ji,jj) = ( ht_0 (ji,jj ) + ht_0 (ji+1,jj ) & 1211 & + ht_0 (ji,jj+1) + ht_0 (ji+1,jj+1) ) & 1212 & / MAX( ssmask(ji,jj ) + ssmask(ji+1,jj ) & 1213 & + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp ) 1214 END DO 1215 END DO 1216 ENDIF 1217 ! 1218 DO jj = 1, jpjm1 1219 zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 1220 END DO 1221 ! 1222 DO jk = 1, jpkm1 1223 DO jj = 1, jpjm1 1224 zhf(:,jj) = zhf(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 1225 END DO 1226 END DO 1227 CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp ) 1228 ! JC: TBC. hf should be greater than 0 1229 DO jj = 1, jpj 1230 DO ji = 1, jpi 1231 IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj) 1232 END DO 1233 END DO 1234 zwz(:,:) = ff_f(:,:) * zwz(:,:) 1235 END SELECT 1236 1237 END SUBROUTINE dyn_cor_2d_init 1238 1239 1240 1241 SUBROUTINE dyn_cor_2d( ht_n, hu_n, hv_n, un_b, vn_b, zhU, zhV, zu_trd, zv_trd ) 1242 !!--------------------------------------------------------------------- 1243 !! *** ROUTINE dyn_cor_2d *** 1244 !! 1245 !! ** Purpose : Compute u and v coriolis trends 1246 !!---------------------------------------------------------------------- 1247 INTEGER :: ji ,jj ! dummy loop indices 1248 REAL(wp) :: zx1, zx2, zy1, zy2 ! - - 1249 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: ht_n, hu_n, hv_n, un_b, vn_b, zhU, zhV 1250 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: zu_trd, zv_trd 1251 !!---------------------------------------------------------------------- 1252 SELECT CASE( nvor_scheme ) 1253 CASE( np_ENT ) ! enstrophy conserving scheme (f-point) 1254 DO jj = 2, jpjm1 1255 DO ji = 2, jpim1 1256 zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * r1_hu_n(ji,jj) & 1257 & * ( e1e2t(ji+1,jj)*ht_n(ji+1,jj)*ff_t(ji+1,jj) * ( vn_b(ji+1,jj) + vn_b(ji+1,jj-1) ) & 1258 & + e1e2t(ji ,jj)*ht_n(ji ,jj)*ff_t(ji ,jj) * ( vn_b(ji ,jj) + vn_b(ji ,jj-1) ) ) 1259 ! 1260 zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * r1_hv_n(ji,jj) & 1261 & * ( e1e2t(ji,jj+1)*ht_n(ji,jj+1)*ff_t(ji,jj+1) * ( un_b(ji,jj+1) + un_b(ji-1,jj+1) ) & 1262 & + e1e2t(ji,jj )*ht_n(ji,jj )*ff_t(ji,jj ) * ( un_b(ji,jj ) + un_b(ji-1,jj ) ) ) 1263 END DO 1264 END DO 1265 ! 1266 CASE( np_ENE , np_MIX ) ! energy conserving scheme (t-point) ENE or MIX 1267 DO jj = 2, jpjm1 1268 DO ji = fs_2, fs_jpim1 ! vector opt. 1269 zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj) 1270 zy2 = ( zhV(ji,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj) 1271 zx1 = ( zhU(ji-1,jj) + zhU(ji-1,jj+1) ) * r1_e2v(ji,jj) 1272 zx2 = ( zhU(ji ,jj) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj) 1273 ! energy conserving formulation for planetary vorticity term 1274 zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 1275 zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 1276 END DO 1277 END DO 1278 ! 1279 CASE( np_ENS ) ! enstrophy conserving scheme (f-point) 1280 DO jj = 2, jpjm1 1281 DO ji = fs_2, fs_jpim1 ! vector opt. 1282 zy1 = r1_8 * ( zhV(ji ,jj-1) + zhV(ji+1,jj-1) & 1283 & + zhV(ji ,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj) 1284 zx1 = - r1_8 * ( zhU(ji-1,jj ) + zhU(ji-1,jj+1) & 1285 & + zhU(ji ,jj ) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj) 1286 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 1287 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 1288 END DO 1289 END DO 1290 ! 1291 CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f) 1292 DO jj = 2, jpjm1 1293 DO ji = fs_2, fs_jpim1 ! vector opt. 1294 zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zhV(ji ,jj ) & 1295 & + ftnw(ji+1,jj) * zhV(ji+1,jj ) & 1296 & + ftse(ji,jj ) * zhV(ji ,jj-1) & 1297 & + ftsw(ji+1,jj) * zhV(ji+1,jj-1) ) 1298 zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zhU(ji-1,jj+1) & 1299 & + ftse(ji,jj+1) * zhU(ji ,jj+1) & 1300 & + ftnw(ji,jj ) * zhU(ji-1,jj ) & 1301 & + ftne(ji,jj ) * zhU(ji ,jj ) ) 1302 END DO 1303 END DO 1304 ! 1305 END SELECT 1306 ! 1307 END SUBROUTINE dyn_cor_2D 1308 1309 1310 SUBROUTINE wad_tmsk( pssh, ptmsk ) 1311 !!---------------------------------------------------------------------- 1312 !! *** ROUTINE wad_lmt *** 1313 !! 1314 !! ** Purpose : set wetting & drying mask at tracer points 1315 !! for the current barotropic sub-step 1316 !! 1317 !! ** Method : ??? 1318 !! 1319 !! ** Action : ptmsk : wetting & drying t-mask 1320 !!---------------------------------------------------------------------- 1321 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pssh ! 1322 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: ptmsk ! 1323 ! 1324 INTEGER :: ji, jj ! dummy loop indices 1325 !!---------------------------------------------------------------------- 1326 ! 1327 IF( ln_wd_dl_rmp ) THEN 1328 DO jj = 1, jpj 1329 DO ji = 1, jpi 1330 IF ( pssh(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN 1331 ! IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN 1332 ptmsk(ji,jj) = 1._wp 1333 ELSEIF( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN 1334 ptmsk(ji,jj) = TANH( 50._wp*( ( pssh(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )*r_rn_wdmin1) ) 1335 ELSE 1336 ptmsk(ji,jj) = 0._wp 1337 ENDIF 1338 END DO 1339 END DO 1340 ELSE 1341 DO jj = 1, jpj 1342 DO ji = 1, jpi 1343 IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN ; ptmsk(ji,jj) = 1._wp 1344 ELSE ; ptmsk(ji,jj) = 0._wp 1345 ENDIF 1346 END DO 1347 END DO 1348 ENDIF 1349 ! 1350 END SUBROUTINE wad_tmsk 1351 1352 1353 SUBROUTINE wad_Umsk( pTmsk, phU, phV, pu, pv, pUmsk, pVmsk ) 1354 !!---------------------------------------------------------------------- 1355 !! *** ROUTINE wad_lmt *** 1356 !! 1357 !! ** Purpose : set wetting & drying mask at tracer points 1358 !! for the current barotropic sub-step 1359 !! 1360 !! ** Method : ??? 1361 !! 1362 !! ** Action : ptmsk : wetting & drying t-mask 1363 !!---------------------------------------------------------------------- 1364 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pTmsk ! W & D t-mask 1365 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phU, phV, pu, pv ! ocean velocities and transports 1366 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pUmsk, pVmsk ! W & D u- and v-mask 1367 ! 1368 INTEGER :: ji, jj ! dummy loop indices 1369 !!---------------------------------------------------------------------- 1370 ! 1371 DO jj = 1, jpj 1372 DO ji = 1, jpim1 ! not jpi-column 1373 IF ( phU(ji,jj) > 0._wp ) THEN ; pUmsk(ji,jj) = pTmsk(ji ,jj) 1374 ELSE ; pUmsk(ji,jj) = pTmsk(ji+1,jj) 1375 ENDIF 1376 phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj) 1377 pu (ji,jj) = pUmsk(ji,jj)*pu (ji,jj) 1378 END DO 1379 END DO 1380 ! 1381 DO jj = 1, jpjm1 ! not jpj-row 1382 DO ji = 1, jpi 1383 IF ( phV(ji,jj) > 0._wp ) THEN ; pVmsk(ji,jj) = pTmsk(ji,jj ) 1384 ELSE ; pVmsk(ji,jj) = pTmsk(ji,jj+1) 1385 ENDIF 1386 phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj) 1387 pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj) 1388 END DO 1389 END DO 1390 ! 1391 END SUBROUTINE wad_Umsk 1392 1393 1394 SUBROUTINE wad_spg( sshn, zcpx, zcpy ) 1395 !!--------------------------------------------------------------------- 1396 !! *** ROUTINE wad_sp *** 1397 !! 1398 !! ** Purpose : 1399 !!---------------------------------------------------------------------- 1400 INTEGER :: ji ,jj ! dummy loop indices 1401 LOGICAL :: ll_tmp1, ll_tmp2 1402 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: sshn 1403 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy 1404 !!---------------------------------------------------------------------- 1405 DO jj = 2, jpjm1 1406 DO ji = 2, jpim1 1407 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 1408 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 1409 & MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji+1,jj) + ht_0(ji+1,jj) ) & 1410 & > rn_wdmin1 + rn_wdmin2 1411 ll_tmp2 = ( ABS( sshn(ji+1,jj) - sshn(ji ,jj)) > 1.E-12 ).AND.( & 1412 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 1413 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 1414 IF(ll_tmp1) THEN 1415 zcpx(ji,jj) = 1.0_wp 1416 ELSEIF(ll_tmp2) THEN 1417 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 1418 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 1419 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 1420 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 1421 ELSE 1422 zcpx(ji,jj) = 0._wp 1423 ENDIF 1424 ! 1425 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 1426 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 1427 & MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji,jj+1) + ht_0(ji,jj+1) ) & 1428 & > rn_wdmin1 + rn_wdmin2 1429 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1)) > 1.E-12 ).AND.( & 1430 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 1431 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 1432 1433 IF(ll_tmp1) THEN 1434 zcpy(ji,jj) = 1.0_wp 1435 ELSE IF(ll_tmp2) THEN 1436 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 1437 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 1438 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 1439 zcpy(ji,jj) = MAX( 0._wp , MIN( zcpy(ji,jj) , 1.0_wp ) ) 1440 ELSE 1441 zcpy(ji,jj) = 0._wp 1442 ENDIF 1443 END DO 1444 END DO 1445 1446 END SUBROUTINE wad_spg 1447 1448 1449 1450 SUBROUTINE drg_init( pu_RHSi, pv_RHSi, pCdU_u, pCdU_v ) 1451 !!---------------------------------------------------------------------- 1452 !! *** ROUTINE drg_init *** 1453 !! 1454 !! ** Purpose : - add the baroclinic top/bottom drag contribution to 1455 !! the baroclinic part of the barotropic RHS 1456 !! - compute the barotropic drag coefficients 1457 !! 1458 !! ** Method : computation done over the INNER domain only 1459 !!---------------------------------------------------------------------- 1460 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pu_RHSi, pv_RHSi ! baroclinic part of the barotropic RHS 1461 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: pCdU_u , pCdU_v ! barotropic drag coefficients 1462 ! 1463 INTEGER :: ji, jj ! dummy loop indices 1464 INTEGER :: ikbu, ikbv, iktu, iktv 1465 REAL(wp) :: zztmp 1466 REAL(wp), DIMENSION(jpi,jpj) :: zu_i, zv_i 1467 !!---------------------------------------------------------------------- 1468 ! 1469 ! !== Set the barotropic drag coef. ==! 1470 ! 1471 IF( ln_isfcav ) THEN ! top+bottom friction (ocean cavities) 1472 1473 DO jj = 2, jpjm1 1474 DO ji = 2, jpim1 ! INNER domain 1475 pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) 1476 pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) 1477 END DO 1478 END DO 1479 ELSE ! bottom friction only 1480 DO jj = 2, jpjm1 1481 DO ji = 2, jpim1 ! INNER domain 1482 pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) 1483 pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) 1484 END DO 1485 END DO 1486 ENDIF 1487 ! 1488 ! !== BOTTOM stress contribution from baroclinic velocities ==! 1489 ! 1490 IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW bottom baroclinic velocities 1491 1492 DO jj = 2, jpjm1 1493 DO ji = 2, jpim1 ! INNER domain 1494 ikbu = mbku(ji,jj) 1495 ikbv = mbkv(ji,jj) 1496 zu_i(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) 1497 zv_i(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj) 1498 END DO 1499 END DO 1500 ELSE ! CENTRED integration: use BEFORE bottom baroclinic velocities 1501 1502 DO jj = 2, jpjm1 1503 DO ji = 2, jpim1 ! INNER domain 1504 ikbu = mbku(ji,jj) 1505 ikbv = mbkv(ji,jj) 1506 zu_i(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) 1507 zv_i(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj) 1508 END DO 1509 END DO 1510 ENDIF 1511 ! 1512 IF( ln_wd_il ) THEN ! W/D : use the "clipped" bottom friction !!gm explain WHY, please ! 1513 zztmp = -1._wp / rdtbt 1514 DO jj = 2, jpjm1 1515 DO ji = 2, jpim1 ! INNER domain 1516 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( & 1517 & r1_hu_n(ji,jj) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp ) 1518 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) * wdrampv(ji,jj) * MAX( & 1519 & r1_hv_n(ji,jj) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp ) 1520 END DO 1521 END DO 1522 ELSE ! use "unclipped" drag (even if explicit friction is used in 3D calculation) 1523 1524 DO jj = 2, jpjm1 1525 DO ji = 2, jpim1 ! INNER domain 1526 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu_n(ji,jj) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zu_i(ji,jj) 1527 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv_n(ji,jj) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zv_i(ji,jj) 1528 END DO 1529 END DO 1530 END IF 1531 ! 1532 ! !== TOP stress contribution from baroclinic velocities ==! (no W/D case) 1533 ! 1534 IF( ln_isfcav ) THEN 1535 ! 1536 IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW top baroclinic velocity 1537 1538 DO jj = 2, jpjm1 1539 DO ji = 2, jpim1 ! INNER domain 1540 iktu = miku(ji,jj) 1541 iktv = mikv(ji,jj) 1542 zu_i(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) 1543 zv_i(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj) 1544 END DO 1545 END DO 1546 ELSE ! CENTRED integration: use BEFORE top baroclinic velocity 1547 1548 DO jj = 2, jpjm1 1549 DO ji = 2, jpim1 ! INNER domain 1550 iktu = miku(ji,jj) 1551 iktv = mikv(ji,jj) 1552 zu_i(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) 1553 zv_i(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj) 1554 END DO 1555 END DO 1556 ENDIF 1557 ! 1558 ! ! use "unclipped" top drag (even if explicit friction is used in 3D calculation) 1559 1560 DO jj = 2, jpjm1 1561 DO ji = 2, jpim1 ! INNER domain 1562 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu_n(ji,jj) * r1_2*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zu_i(ji,jj) 1563 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv_n(ji,jj) * r1_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zv_i(ji,jj) 1564 END DO 1565 END DO 1566 ! 1567 ENDIF 1568 ! 1569 END SUBROUTINE drg_init 1570 1571 SUBROUTINE ts_bck_interp( jn, ll_init, & ! <== in 1572 & za0, za1, za2, za3 ) ! ==> out 1573 !!---------------------------------------------------------------------- 1574 INTEGER ,INTENT(in ) :: jn ! index of sub time step 1575 LOGICAL ,INTENT(in ) :: ll_init ! 1576 REAL(wp),INTENT( out) :: za0, za1, za2, za3 ! Half-step back interpolation coefficient 1577 ! 1578 REAL(wp) :: zepsilon, zgamma ! - - 1579 !!---------------------------------------------------------------------- 1580 ! ! set Half-step back interpolation coefficient 1581 IF ( jn==1 .AND. ll_init ) THEN !* Forward-backward 1582 za0 = 1._wp 1583 za1 = 0._wp 1584 za2 = 0._wp 1585 za3 = 0._wp 1586 ELSEIF( jn==2 .AND. ll_init ) THEN !* AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12 1587 za0 = 1.0833333333333_wp ! za0 = 1-gam-eps 1588 za1 =-0.1666666666666_wp ! za1 = gam 1589 za2 = 0.0833333333333_wp ! za2 = eps 1590 za3 = 0._wp 1591 ELSE !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880 1592 IF( rn_bt_alpha == 0._wp ) THEN ! Time diffusion 1593 za0 = 0.614_wp ! za0 = 1/2 + gam + 2*eps 1594 za1 = 0.285_wp ! za1 = 1/2 - 2*gam - 3*eps 1595 za2 = 0.088_wp ! za2 = gam 1596 za3 = 0.013_wp ! za3 = eps 1597 ELSE ! no time diffusion 1598 zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha 1599 zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha 1600 za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon 1601 za1 = 1._wp - za0 - zgamma - zepsilon 1602 za2 = zgamma 1603 za3 = zepsilon 1604 ENDIF 1605 ENDIF 1606 END SUBROUTINE ts_bck_interp 1607 1608 1581 1609 !!====================================================================== 1582 1610 END MODULE dynspg_ts
Note: See TracChangeset
for help on using the changeset viewer.
