Changeset 13231 for NEMO/branches/2020
- Timestamp:
- 2020-07-02T18:03:12+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/NST/agrif_oce_sponge.F90
r13230 r13231 32 32 33 33 !! * Substitutions 34 # include "do_loop_substitute.h90" 35 !!---------------------------------------------------------------------- 36 !! NEMO/NST 4.0 , NEMO Consortium (2018) 37 !! $Id$ 38 !! Software governed by the CeCILL license (see ./LICENSE) 39 !!---------------------------------------------------------------------- 40 CONTAINS 41 42 SUBROUTINE Agrif_Sponge_Tra 43 !!---------------------------------------------------------------------- 44 !! *** ROUTINE Agrif_Sponge_Tra *** 45 !!---------------------------------------------------------------------- 46 REAL(wp) :: zcoef ! local scalar 47 48 !!---------------------------------------------------------------------- 49 ! 50 #if defined SPONGE 51 !! Assume persistence: 52 zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot()) 53 54 CALL Agrif_Sponge 55 Agrif_SpecialValue = 0._wp 56 Agrif_UseSpecialValue = .TRUE. 57 tabspongedone_tsn = .FALSE. 58 ! 59 CALL Agrif_Bc_Variable( tsn_sponge_id, calledweight=zcoef, procname=interptsn_sponge ) 60 ! 61 Agrif_UseSpecialValue = .FALSE. 62 #endif 63 ! 64 CALL iom_put( 'agrif_spu', fspu(:,:)) 65 CALL iom_put( 'agrif_spv', fspv(:,:)) 66 ! 67 END SUBROUTINE Agrif_Sponge_Tra 68 69 70 SUBROUTINE Agrif_Sponge_dyn 71 !!---------------------------------------------------------------------- 72 !! *** ROUTINE Agrif_Sponge_dyn *** 73 !!---------------------------------------------------------------------- 74 REAL(wp) :: zcoef ! local scalar 75 !!---------------------------------------------------------------------- 76 ! 77 #if defined SPONGE 78 zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot()) 79 80 Agrif_SpecialValue=0. 81 Agrif_UseSpecialValue = ln_spc_dyn 82 use_sign_north = .TRUE. 83 sign_north = -1. 84 ! 85 tabspongedone_u = .FALSE. 86 tabspongedone_v = .FALSE. 87 CALL Agrif_Bc_Variable( un_sponge_id, calledweight=zcoef, procname=interpun_sponge ) 88 ! 89 tabspongedone_u = .FALSE. 90 tabspongedone_v = .FALSE. 91 CALL Agrif_Bc_Variable( vn_sponge_id, calledweight=zcoef, procname=interpvn_sponge ) 92 ! 93 Agrif_UseSpecialValue = .FALSE. 94 use_sign_north = .FALSE. 95 #endif 96 ! 97 CALL iom_put( 'agrif_spt', fspt(:,:)) 98 CALL iom_put( 'agrif_spf', fspf(:,:)) 99 ! 100 END SUBROUTINE Agrif_Sponge_dyn 101 102 103 SUBROUTINE Agrif_Sponge 104 !!---------------------------------------------------------------------- 105 !! *** ROUTINE Agrif_Sponge *** 106 !!---------------------------------------------------------------------- 107 INTEGER :: ji, jj, ind1, ind2 108 INTEGER :: ispongearea, jspongearea 109 REAL(wp) :: z1_ispongearea, z1_jspongearea 110 REAL(wp), DIMENSION(jpi,jpj) :: ztabramp 111 #if defined key_vertical 112 REAL(wp), DIMENSION(jpi,jpj) :: ztabrampu 113 REAL(wp), DIMENSION(jpi,jpj) :: ztabrampv 114 #endif 115 REAL(wp), DIMENSION(jpjmax) :: zmskwest, zmskeast 116 REAL(wp), DIMENSION(jpimax) :: zmsknorth, zmsksouth 117 !!---------------------------------------------------------------------- 118 ! 119 ! Sponge 1d example with: 120 ! iraf = 3 ; nbghost = 3 ; nn_sponge_len = 2 121 ! 122 !coarse : U T U T U T U 123 !| | | | | 124 !fine : t u t u t u t u t u t u t u t u t u t u t 125 !sponge val:0 0 0 1 5/6 4/6 3/6 2/6 1/6 0 0 126 ! | ghost | <-- sponge area -- > | 127 ! | points | | 128 ! |--> dynamical interface 129 130 #if defined SPONGE || defined SPONGE_TOP 131 IF (( .NOT. spongedoneT ).OR.( .NOT. spongedoneU )) THEN 132 ! 133 ! Retrieve masks at open boundaries: 134 135 ! --- West --- ! 136 IF( lk_west) THEN 137 ztabramp(:,:) = 0._wp 138 ind1 = nn_hls + 1 + nbghostcells ! halo + land + nbghostcells 139 DO ji = mi0(ind1), mi1(ind1) 140 ztabramp(ji,:) = ssumask(ji,:) 141 END DO 142 zmskwest( 1:jpj ) = MAXVAL(ztabramp(:,:), dim=1) 143 zmskwest(jpj+1:jpjmax) = 0._wp 144 ENDIF 145 146 ! --- East --- ! 147 IF( lk_east ) THEN 148 ztabramp(:,:) = 0._wp 149 ind1 = jpiglo - ( nn_hls + nbghostcells + 1) ! halo + land + nbghostcells 150 DO ji = mi0(ind1), mi1(ind1) 151 ztabramp(ji,:) = ssumask(ji,:) 152 END DO 153 zmskeast( 1:jpj ) = MAXVAL(ztabramp(:,:), dim=1) 154 zmskeast(jpj+1:jpjmax) = 0._wp 155 ENDIF 156 157 ! --- South --- ! 158 IF( lk_south ) THEN 159 ztabramp(:,:) = 0._wp 160 ind1 = nn_hls + 1 + nbghostcells ! halo + land + nbghostcells 161 DO jj = mj0(ind1), mj1(ind1) 162 ztabramp(:,jj) = ssvmask(:,jj) 163 END DO 164 zmsksouth( 1:jpi ) = MAXVAL(ztabramp(:,:), dim=2) 165 zmsksouth(jpi+1:jpimax) = 0._wp 166 ENDIF 167 168 ! --- North --- ! 169 IF( lk_north) THEN 170 ztabramp(:,:) = 0._wp 171 ind1 = jpjglo - ( nn_hls + nbghostcells + 1) ! halo + land + nbghostcells 172 DO jj = mj0(ind1), mj1(ind1) 173 ztabramp(:,jj) = ssvmask(:,jj) 174 END DO 175 zmsknorth( 1:jpi ) = MAXVAL(ztabramp(:,:), dim=2) 176 zmsknorth(jpi+1:jpimax) = 0._wp 177 ENDIF 178 179 ! JC: SPONGE MASKING TO BE SORTED OUT: 180 zmskwest(:) = 1._wp 181 zmskeast(:) = 1._wp 182 zmsksouth(:) = 1._wp 183 zmsknorth(:) = 1._wp 184 #if defined key_mpp_mpi 185 ! CALL mpp_max( 'AGRIF_sponge', zmskwest(:) , jpjmax ) 186 ! CALL mpp_max( 'AGRIF_Sponge', zmskeast(:) , jpjmax ) 187 ! CALL mpp_max( 'AGRIF_Sponge', zmsksouth(:), jpimax ) 188 ! CALL mpp_max( 'AGRIF_Sponge', zmsknorth(:), jpimax ) 189 #endif 190 191 ! Define ramp from boundaries towards domain interior at T-points 192 ! Store it in ztabramp 193 194 ispongearea = nn_sponge_len * Agrif_irhox() 195 z1_ispongearea = 1._wp / REAL( ispongearea ) 196 jspongearea = nn_sponge_len * Agrif_irhoy() 197 z1_jspongearea = 1._wp / REAL( jspongearea ) 198 199 ztabramp(:,:) = 0._wp 200 201 ! Trick to remove sponge in 2DV domains: 202 IF ( nbcellsx <= 3 ) ispongearea = -1 203 IF ( nbcellsy <= 3 ) jspongearea = -1 204 205 ! --- West --- ! 206 IF(lk_west) THEN 207 ind1 = nn_hls + 1 + nbghostcells ! halo + land + nbghostcells 208 ind2 = nn_hls + 1 + nbghostcells + ispongearea 209 DO ji = mi0(ind1), mi1(ind2) 210 DO jj = 1, jpj 211 ztabramp(ji,jj) = REAL( ind2 - mig(ji) ) * z1_ispongearea * zmskwest(jj) 212 END DO 213 END DO 214 ! ghost cells: 215 ind1 = 1 216 ind2 = nn_hls + 1 + nbghostcells ! halo + land + nbghostcells 217 DO ji = mi0(ind1), mi1(ind2) 218 DO jj = 1, jpj 219 ztabramp(ji,jj) = zmskwest(jj) 220 END DO 221 END DO 222 ENDIF 223 224 ! --- East --- ! 225 IF(lk_east) THEN 226 ind1 = jpiglo - ( nn_hls + nbghostcells ) - ispongearea 227 ind2 = jpiglo - ( nn_hls + nbghostcells ) ! halo + land + nbghostcells - 1 228 DO ji = mi0(ind1), mi1(ind2) 229 DO jj = 1, jpj 230 ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( mig(ji) - ind1 ) * z1_ispongearea) * zmskeast(jj) 231 END DO 232 END DO 233 ! ghost cells: 234 ind1 = jpiglo - ( nn_hls + nbghostcells ) ! halo + land + nbghostcells - 1 235 ind2 = jpiglo 236 DO ji = mi0(ind1), mi1(ind2) 237 DO jj = 1, jpj 238 ztabramp(ji,jj) = zmskeast(jj) 239 END DO 240 END DO 241 ENDIF 242 243 ! --- South --- ! 244 IF( lk_south ) THEN 245 ind1 = nn_hls + 1 + nbghostcells ! halo + land + nbghostcells 246 ind2 = nn_hls + 1 + nbghostcells + jspongearea 247 DO jj = mj0(ind1), mj1(ind2) 248 DO ji = 1, jpi 249 ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( ind2 - mjg(jj) ) * z1_jspongearea) * zmsksouth(ji) 250 END DO 251 END DO 252 ! ghost cells: 253 ind1 = 1 254 ind2 = nn_hls + 1 + nbghostcells ! halo + land + nbghostcells 255 DO jj = mj0(ind1), mj1(ind2) 256 DO ji = 1, jpi 257 ztabramp(ji,jj) = zmsksouth(ji) 258 END DO 259 END DO 260 ENDIF 261 262 ! --- North --- ! 263 IF( lk_north ) THEN 264 ind1 = jpjglo - ( nn_hls + nbghostcells ) - jspongearea 265 ind2 = jpjglo - ( nn_hls + nbghostcells ) ! halo + land + nbghostcells - 1 266 DO jj = mj0(ind1), mj1(ind2) 267 DO ji = 1, jpi 268 ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( mjg(jj) - ind1 ) * z1_jspongearea) * zmsknorth(ji) 269 END DO 270 END DO 271 ! ghost cells: 272 ind1 = jpjglo - ( nn_hls + nbghostcells ) ! halo + land + nbghostcells - 1 273 ind2 = jpjglo 274 DO jj = mj0(ind1), mj1(ind2) 275 DO ji = 1, jpi 276 ztabramp(ji,jj) = zmsknorth(ji) 277 END DO 278 END DO 279 ENDIF 280 ! 281 ENDIF 282 283 ! Tracers 284 IF( .NOT. spongedoneT ) THEN 285 fspu(:,:) = 0._wp 286 fspv(:,:) = 0._wp 287 DO_2D_00_00 288 fspu(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji+1,jj ) ) * ssumask(ji,jj) 289 fspv(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji ,jj+1) ) * ssvmask(ji,jj) 290 END_2D 291 ENDIF 292 293 ! Dynamics 294 IF( .NOT. spongedoneU ) THEN 295 fspt(:,:) = 0._wp 296 fspf(:,:) = 0._wp 297 DO_2D_00_00 298 fspt(ji,jj) = ztabramp(ji,jj) * ssmask(ji,jj) 299 fspf(ji,jj) = 0.25_wp * ( ztabramp(ji ,jj ) + ztabramp(ji ,jj+1) & 300 & +ztabramp(ji+1,jj+1) + ztabramp(ji+1,jj ) ) & 301 & * ssvmask(ji,jj) * ssvmask(ji,jj+1) 302 END_2D 303 ENDIF 304 305 IF( .NOT. spongedoneT .AND. .NOT. spongedoneU ) THEN 306 CALL lbc_lnk_multi( 'agrif_Sponge', fspu, 'U', 1., fspv, 'V', 1., fspt, 'T', 1., fspf, 'F', 1. ) 307 spongedoneT = .TRUE. 308 spongedoneU = .TRUE. 309 ENDIF 310 IF( .NOT. spongedoneT ) THEN 311 CALL lbc_lnk_multi( 'agrif_Sponge', fspu, 'U', 1., fspv, 'V', 1. ) 312 spongedoneT = .TRUE. 313 ENDIF 314 IF( .NOT. spongedoneT .AND. .NOT. spongedoneU ) THEN 315 CALL lbc_lnk_multi( 'agrif_Sponge', fspt, 'T', 1., fspf, 'F', 1. ) 316 spongedoneU = .TRUE. 317 ENDIF 318 319 #if defined key_vertical 320 ! Remove vertical interpolation where not needed: 321 DO_2D_00_00 322 IF ((fspu(ji-1,jj)==0._wp).AND.(fspu(ji,jj)==0._wp).AND. & 323 & (fspv(ji,jj-1)==0._wp).AND.(fspv(ji,jj)==0._wp)) mbkt_parent(ji,jj) = 0 324 ! 325 IF ((fspt(ji+1,jj)==0._wp).AND.(fspt(ji,jj)==0._wp).AND. & 326 & (fspf(ji,jj-1)==0._wp).AND.(fspf(ji,jj)==0._wp)) mbku_parent(ji,jj) = 0 327 ! 328 IF ((fspt(ji,jj+1)==0._wp).AND.(fspt(ji,jj)==0._wp).AND. & 329 & (fspf(ji-1,jj)==0._wp).AND.(fspf(ji,jj)==0._wp)) mbkv_parent(ji,jj) = 0 330 ! 331 IF ( ssmask(ji,jj) == 0._wp) mbkt_parent(ji,jj) = 0 332 IF (ssumask(ji,jj) == 0._wp) mbku_parent(ji,jj) = 0 333 IF (ssvmask(ji,jj) == 0._wp) mbkv_parent(ji,jj) = 0 334 END_2D 335 ! 336 ztabramp (:,:) = REAL( mbkt_parent (:,:), wp ) 337 ztabrampu(:,:) = REAL( mbku_parentu(:,:), wp ) 338 ztabrampv(:,:) = REAL( mbkv_parentv(:,:), wp ) 339 CALL lbc_lnk_multi( 'Agrif_Sponge', ztabramp, 'T', 1., ztabrampu, 'U', 1., ztabrampv, 'V', 1. ) 340 mbkt_parent(:,:) = NINT( ztabramp (:,:) ) 341 mbku_parent(:,:) = NINT( ztabrampu(:,:) ) 342 mbkv_parent(:,:) = NINT( ztabrampv(:,:) ) 343 #endif 344 ! 345 #endif 346 ! 347 END SUBROUTINE Agrif_Sponge 348 349 SUBROUTINE interptsn_sponge( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 350 !!---------------------------------------------------------------------- 351 !! *** ROUTINE interptsn_sponge *** 352 !!---------------------------------------------------------------------- 353 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2 354 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres 355 LOGICAL , INTENT(in ) :: before 356 ! 357 INTEGER :: ji, jj, jk, jn ! dummy loop indices 358 INTEGER :: iku, ikv 359 REAL(wp) :: ztsa, zabe1, zabe2, zbtr, zhtot 360 REAL(wp), DIMENSION(i1:i2,j1:j2,jpk) :: ztu, ztv 361 REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tsbdiff 362 ! vertical interpolation: 363 REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tabres_child 364 REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin 365 REAL(wp), DIMENSION(k1:k2) :: h_in 366 REAL(wp), DIMENSION(1:jpk) :: h_out 367 INTEGER :: N_in, N_out 368 !!---------------------------------------------------------------------- 369 ! 370 IF( before ) THEN 371 DO jn = 1, jpts 372 DO jk=k1,k2 373 DO jj=j1,j2 374 DO ji=i1,i2 375 tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kbb_a) 376 END DO 377 END DO 378 END DO 379 END DO 380 381 # if defined key_vertical 382 ! Interpolate thicknesses 383 ! Warning: these are masked, hence extrapolated prior interpolation. 384 DO jk=k1,k2 385 DO jj=j1,j2 386 DO ji=i1,i2 387 tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kbb_a) 388 END DO 389 END DO 390 END DO 391 392 ! Extrapolate thicknesses in partial bottom cells: 393 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 394 IF (ln_zps) THEN 395 DO jj=j1,j2 396 DO ji=i1,i2 397 jk = mbkt(ji,jj) 398 tabres(ji,jj,jk,jpts+1) = 0._wp 399 END DO 400 END DO 401 END IF 402 403 ! Save ssh at last level: 404 IF (.NOT.ln_linssh) THEN 405 tabres(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kbb_a)*tmask(i1:i2,j1:j2,1) 406 ELSE 407 tabres(i1:i2,j1:j2,k2,jpts+1) = 0._wp 408 END IF 409 # endif 410 411 ELSE 412 ! 413 # if defined key_vertical 414 415 IF (ln_linssh) tabres(i1:i2,j1:j2,k2,n2) = 0._wp 416 417 DO jj=j1,j2 418 DO ji=i1,i2 419 tabres_child(ji,jj,:,:) = 0._wp 420 N_in = mbkt_parent(ji,jj) 421 zhtot = 0._wp 422 DO jk=1,N_in !k2 = jpk of parent grid 423 IF (jk==N_in) THEN 424 h_in(jk) = ht0_parent(ji,jj) + tabres(ji,jj,k2,n2) - zhtot 425 ELSE 426 h_in(jk) = tabres(ji,jj,jk,n2) 427 ENDIF 428 zhtot = zhtot + h_in(jk) 429 tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1) 430 END DO 431 N_out = 0 432 DO jk=1,jpk ! jpk of child grid 433 IF (tmask(ji,jj,jk) == 0) EXIT 434 N_out = N_out + 1 435 h_out(jk) = e3t(ji,jj,jk,Kbb_a) !Child grid scale factors. Could multiply by e1e2t here instead of division above 436 END DO 437 438 ! Account for small differences in free-surface 439 IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN 440 h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) ) 441 ELSE 442 h_in(1) = h_in(1) - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) ) 443 ENDIF 444 IF (N_in*N_out > 0) THEN 445 CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts) 446 ENDIF 447 END DO 448 END DO 449 # endif 450 451 DO jj=j1,j2 452 DO ji=i1,i2 453 DO jk=1,jpkm1 454 # if defined key_vertical 455 tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres_child(ji,jj,jk,1:jpts)) * tmask(ji,jj,jk) 456 # else 457 tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres(ji,jj,jk,1:jpts))*tmask(ji,jj,jk) 458 # endif 459 END DO 460 END DO 461 END DO 462 463 DO jn = 1, jpts 464 DO jk = 1, jpkm1 465 ztu(i1:i2,j1:j2,jk) = 0._wp 466 DO jj = j1,j2 467 DO ji = i1,i2-1 468 zabe1 = rn_sponge_tra * r1_Dt * fspu(ji,jj) * umask(ji,jj,jk) * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) 469 ztu(ji,jj,jk) = zabe1 * ( tsbdiff(ji+1,jj ,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 470 END DO 471 END DO 472 ztv(i1:i2,j1:j2,jk) = 0._wp 473 DO ji = i1,i2 474 DO jj = j1,j2-1 475 zabe2 = rn_sponge_tra * r1_Dt * fspv(ji,jj) * vmask(ji,jj,jk) * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a) 476 ztv(ji,jj,jk) = zabe2 * ( tsbdiff(ji ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 477 END DO 478 END DO 479 ! 480 IF( ln_zps ) THEN ! set gradient at partial step level 481 DO jj = j1,j2 482 DO ji = i1,i2 483 ! last level 484 iku = mbku(ji,jj) 485 ikv = mbkv(ji,jj) 486 IF( iku == jk ) ztu(ji,jj,jk) = 0._wp 487 IF( ikv == jk ) ztv(ji,jj,jk) = 0._wp 488 END DO 489 END DO 490 ENDIF 491 END DO 492 ! 493 DO jk = 1, jpkm1 494 DO jj = j1+1,j2-1 495 DO ji = i1+1,i2-1 496 IF (.NOT. tabspongedone_tsn(ji,jj)) THEN 497 zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm_a) 498 ! horizontal diffusive trends 499 ztsa = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) & 500 & - rn_trelax_tra * r1_Dt * fspt(ji,jj) * tsbdiff(ji,jj,jk,jn) 501 ! add it to the general tracer trends 502 ts(ji,jj,jk,jn,Krhs_a) = ts(ji,jj,jk,jn,Krhs_a) + ztsa 503 ENDIF 504 END DO 505 END DO 506 END DO 507 ! 508 END DO 509 ! 510 tabspongedone_tsn(i1+1:i2-1,j1+1:j2-1) = .TRUE. 511 ! 512 ENDIF 513 ! 514 END SUBROUTINE interptsn_sponge 515 516 SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before) 517 !!--------------------------------------------- 518 !! *** ROUTINE interpun_sponge *** 519 !!--------------------------------------------- 520 INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2 521 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres 522 LOGICAL, INTENT(in) :: before 523 524 INTEGER :: ji,jj,jk,jmax 525 INTEGER :: ind1 526 ! sponge parameters 527 REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, zhtot 528 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ubdiff 529 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff 530 ! vertical interpolation: 531 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child 532 REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 533 REAL(wp), DIMENSION(1:jpk) :: h_out 534 INTEGER ::N_in, N_out 535 !!--------------------------------------------- 536 ! 537 IF( before ) THEN 538 DO jk=k1,k2 539 DO jj=j1,j2 540 DO ji=i1,i2 541 tabres(ji,jj,jk,m1) = uu(ji,jj,jk,Kbb_a) 542 # if defined key_vertical 543 tabres(ji,jj,jk,m2) = e3u(ji,jj,jk,Kbb_a)*umask(ji,jj,jk) 544 # endif 545 END DO 546 END DO 547 END DO 548 549 # if defined key_vertical 550 ! Extrapolate thicknesses in partial bottom cells: 551 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 552 IF (ln_zps) THEN 553 DO jj=j1,j2 554 DO ji=i1,i2 555 jk = mbku(ji,jj) 556 tabres(ji,jj,jk,m2) = 0._wp 557 END DO 558 END DO 559 END IF 560 ! Save ssh at last level: 561 tabres(i1:i2,j1:j2,k2,m2) = 0._wp 562 IF (.NOT.ln_linssh) THEN 563 ! This vertical sum below should be replaced by the sea-level at U-points (optimization): 564 DO jk=1,jpk 565 tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) + e3u(i1:i2,j1:j2,jk,Kbb_a) * umask(i1:i2,j1:j2,jk) 566 END DO 567 tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hu_0(i1:i2,j1:j2) 568 END IF 569 # endif 570 571 ELSE 572 573 # if defined key_vertical 574 IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp 575 576 DO jj=j1,j2 577 DO ji=i1,i2 578 tabres_child(ji,jj,:) = 0._wp 579 N_in = mbku_parent(ji,jj) 580 zhtot = 0._wp 581 DO jk=1,N_in 582 IF (jk==N_in) THEN 583 h_in(jk) = hu0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot 584 ELSE 585 h_in(jk) = tabres(ji,jj,jk,m2) 586 ENDIF 587 zhtot = zhtot + h_in(jk) 588 tabin(jk) = tabres(ji,jj,jk,m1) 589 END DO 590 ! 591 N_out = 0 592 DO jk=1,jpk 593 IF (umask(ji,jj,jk) == 0) EXIT 594 N_out = N_out + 1 595 h_out(N_out) = e3u(ji,jj,jk,Kbb_a) 596 END DO 597 598 ! Account for small differences in free-surface 599 IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN 600 h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) ) 601 ELSE 602 h_in(1) = h_in(1) - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) ) 603 ENDIF 604 605 IF (N_in * N_out > 0) THEN 606 CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1) 607 ENDIF 608 END DO 609 END DO 610 611 ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:) 612 #else 613 ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:) 614 #endif 615 ! 616 DO jk = 1, jpkm1 ! Horizontal slab 617 ! ! =============== 618 619 ! ! -------- 620 ! Horizontal divergence ! div 621 ! ! -------- 622 DO jj = j1,j2 623 DO ji = i1+1,i2 ! vector opt. 624 zbtr = rn_sponge_dyn * r1_Dt * fspt(ji,jj) / e3t(ji,jj,jk,Kbb_a) 625 hdivdiff(ji,jj,jk) = ( e2u(ji ,jj)*e3u(ji ,jj,jk,Kbb_a) * ubdiff(ji ,jj,jk) & 626 & -e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb_a) * ubdiff(ji-1,jj,jk) ) * zbtr 627 END DO 628 END DO 629 630 DO jj = j1,j2-1 631 DO ji = i1,i2 ! vector opt. 632 zbtr = rn_sponge_dyn * r1_Dt * fspf(ji,jj) * e3f(ji,jj,jk) 633 rotdiff(ji,jj,jk) = ( -e1u(ji,jj+1) * ubdiff(ji,jj+1,jk) & 634 & +e1u(ji,jj ) * ubdiff(ji,jj ,jk) ) * fmask(ji,jj,jk) * zbtr 635 END DO 636 END DO 637 END DO 638 ! 639 DO jj = j1+1, j2-1 640 DO ji = i1+1, i2-1 ! vector opt. 641 642 IF (.NOT. tabspongedone_u(ji,jj)) THEN 643 DO jk = 1, jpkm1 ! Horizontal slab 644 ze2u = rotdiff (ji,jj,jk) 645 ze1v = hdivdiff(ji,jj,jk) 646 ! horizontal diffusive trends 647 zua = - ( ze2u - rotdiff (ji,jj-1,jk) ) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) ) & 648 & + ( hdivdiff(ji+1,jj,jk) - ze1v ) * r1_e1u(ji,jj) & 649 & - rn_trelax_dyn * r1_Dt * fspu(ji,jj) * ubdiff(ji,jj,jk) 650 651 ! add it to the general momentum trends 652 uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a) + zua 653 END DO 654 ENDIF 655 656 END DO 657 END DO 658 659 tabspongedone_u(i1+1:i2-1,j1+1:j2-1) = .TRUE. 660 661 jmax = j2-1 662 ind1 = jpjglo - ( nn_hls + nbghostcells + 2 ) ! North 663 DO jj = mj0(ind1), mj1(ind1) 664 jmax = MIN(jmax,jj) 665 END DO 666 667 DO jj = j1+1, jmax 668 DO ji = i1+1, i2 ! vector opt. 669 670 IF (.NOT. tabspongedone_v(ji,jj)) THEN 671 DO jk = 1, jpkm1 ! Horizontal slab 672 ze2u = rotdiff (ji,jj,jk) 673 ze1v = hdivdiff(ji,jj,jk) 674 675 ! horizontal diffusive trends 676 zva = + ( ze2u - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) ) & 677 + ( hdivdiff(ji,jj+1,jk) - ze1v ) * r1_e2v(ji,jj) 678 679 ! add it to the general momentum trends 680 vv(ji,jj,jk,Krhs_a) = vv(ji,jj,jk,Krhs_a) + zva 681 END DO 682 ENDIF 683 ! 684 END DO 685 END DO 686 ! 687 tabspongedone_v(i1+1:i2,j1+1:jmax) = .TRUE. 688 ! 689 ENDIF 690 ! 691 END SUBROUTINE interpun_sponge 692 693 SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before,nb,ndir) 694 !!--------------------------------------------- 695 !! *** ROUTINE interpvn_sponge *** 696 !!--------------------------------------------- 697 INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2 698 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres 699 LOGICAL, INTENT(in) :: before 700 INTEGER, INTENT(in) :: nb , ndir 701 ! 702 INTEGER :: ji, jj, jk, imax 703 REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, zhtot 704 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: vbdiff 705 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff 706 ! vertical interpolation: 707 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child 708 REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 709 REAL(wp), DIMENSION(1:jpk) :: h_out 710 INTEGER :: N_in, N_out 711 !!--------------------------------------------- 712 713 IF( before ) THEN 714 DO jk=k1,k2 715 DO jj=j1,j2 34 716 DO ji=i1,i2 35 717 tabres(ji,jj,jk,m1) = vv(ji,jj,jk,Kbb_a)
Note: See TracChangeset
for help on using the changeset viewer.
