Changeset 13230 for NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/NST/agrif_oce_sponge.F90
- Timestamp:
- 2020-07-02T17:50:26+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
r13229 r13230 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 CONTAINS41 42 SUBROUTINE Agrif_Sponge_Tra43 !!----------------------------------------------------------------------44 !! *** ROUTINE Agrif_Sponge_Tra ***45 !!----------------------------------------------------------------------46 REAL(wp) :: zcoef ! local scalar47 48 !!----------------------------------------------------------------------49 !50 #if defined SPONGE51 !! Assume persistence:52 zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())53 54 CALL Agrif_Sponge55 Agrif_SpecialValue = 0._wp56 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 #endif63 !64 CALL iom_put( 'agrif_spu', fspu(:,:))65 CALL iom_put( 'agrif_spv', fspv(:,:))66 !67 END SUBROUTINE Agrif_Sponge_Tra68 69 70 SUBROUTINE Agrif_Sponge_dyn71 !!----------------------------------------------------------------------72 !! *** ROUTINE Agrif_Sponge_dyn ***73 !!----------------------------------------------------------------------74 REAL(wp) :: zcoef ! local scalar75 !!----------------------------------------------------------------------76 !77 #if defined SPONGE78 zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())79 80 Agrif_SpecialValue=0.81 Agrif_UseSpecialValue = ln_spc_dyn82 !83 tabspongedone_u = .FALSE.84 tabspongedone_v = .FALSE.85 CALL Agrif_Bc_Variable( un_sponge_id, calledweight=zcoef, procname=interpun_sponge )86 !87 tabspongedone_u = .FALSE.88 tabspongedone_v = .FALSE.89 CALL Agrif_Bc_Variable( vn_sponge_id, calledweight=zcoef, procname=interpvn_sponge )90 !91 Agrif_UseSpecialValue = .FALSE.92 #endif93 !94 CALL iom_put( 'agrif_spt', fspt(:,:))95 CALL iom_put( 'agrif_spf', fspf(:,:))96 !97 END SUBROUTINE Agrif_Sponge_dyn98 99 100 SUBROUTINE Agrif_Sponge101 !!----------------------------------------------------------------------102 !! *** ROUTINE Agrif_Sponge ***103 !!----------------------------------------------------------------------104 INTEGER :: ji, jj, ind1, ind2105 INTEGER :: ispongearea, jspongearea106 REAL(wp) :: z1_ispongearea, z1_jspongearea107 REAL(wp), DIMENSION(jpi,jpj) :: ztabramp108 #if defined key_vertical109 REAL(wp), DIMENSION(jpi,jpj) :: ztabrampu110 REAL(wp), DIMENSION(jpi,jpj) :: ztabrampv111 #endif112 REAL(wp), DIMENSION(jpjmax) :: zmskwest, zmskeast113 REAL(wp), DIMENSION(jpimax) :: zmsknorth, zmsksouth114 !!----------------------------------------------------------------------115 !116 ! Sponge 1d example with:117 ! iraf = 3 ; nbghost = 3 ; nn_sponge_len = 2118 !119 !coarse : U T U T U T U120 !| | | | |121 !fine : t u t u t u t u t u t u t u t u t u t u t122 !sponge val:0 0 0 1 5/6 4/6 3/6 2/6 1/6 0 0123 ! | ghost | <-- sponge area -- > |124 ! | points | |125 ! |--> dynamical interface126 127 #if defined SPONGE || defined SPONGE_TOP128 IF (( .NOT. spongedoneT ).OR.( .NOT. spongedoneU )) THEN129 !130 ! Retrieve masks at open boundaries:131 132 ! --- West --- !133 IF( lk_west ) THEN134 ztabramp(:,:) = 0._wp135 ind1 = nn_hls + 1 + nbghostcells ! halo + land + nbghostcells136 DO ji = mi0(ind1), mi1(ind1)137 ztabramp(ji,:) = ssumask(ji,:)138 END DO139 !140 zmskwest(1:jpj) = MAXVAL(ztabramp(:,:), dim=1)141 zmskwest(jpj+1:jpjmax) = 0._wp142 ENDIF143 144 ! --- East --- !145 IF( lk_east ) THEN146 ztabramp(:,:) = 0._wp147 ind1 = jpiglo - ( nn_hls + nbghostcells + 1) ! halo + land + nbghostcells148 DO ji = mi0(ind1), mi1(ind1)149 ztabramp(ji,:) = ssumask(ji,:)150 END DO151 !152 zmskeast(1:jpj) = MAXVAL(ztabramp(:,:), dim=1)153 zmskeast(jpj+1:jpjmax) = 0._wp154 ENDIF155 156 ! --- South --- !157 IF( lk_south ) THEN158 ztabramp(:,:) = 0._wp159 ind1 = nn_hls + 1 + nbghostcells ! halo + land + nbghostcells160 DO jj = mj0(ind1), mj1(ind1)161 ztabramp(:,jj) = ssvmask(:,jj)162 END DO163 !164 zmsksouth(1:jpi) = MAXVAL(ztabramp(:,:), dim=2)165 zmsksouth(jpi+1:jpimax) = 0._wp166 ENDIF167 168 ! --- North --- !169 IF( lk_north ) THEN170 ztabramp(:,:) = 0._wp171 ind1 = jpjglo - ( nn_hls + nbghostcells + 1) ! halo + land + nbghostcells172 DO jj = mj0(ind1), mj1(ind1)173 ztabramp(:,jj) = ssvmask(:,jj)174 END DO175 !176 zmsknorth(1:jpi) = MAXVAL(ztabramp(:,:), dim=2)177 zmsknorth(jpi+1:jpimax) = 0._wp178 ENDIF179 180 ! JC: SPONGE MASKING TO BE SORTED OUT:181 zmskwest(:) = 1._wp182 zmskeast(:) = 1._wp183 zmsksouth(:) = 1._wp184 zmsknorth(:) = 1._wp185 #if defined key_mpp_mpi186 ! CALL mpp_max( 'AGRIF_sponge', zmskwest(:) , jpjmax )187 ! CALL mpp_max( 'AGRIF_Sponge', zmskeast(:) , jpjmax )188 ! CALL mpp_max( 'AGRIF_Sponge', zmsksouth(:), jpimax )189 ! CALL mpp_max( 'AGRIF_Sponge', zmsknorth(:), jpimax )190 #endif191 192 ! Define ramp from boundaries towards domain interior at T-points193 ! Store it in ztabramp194 195 ispongearea = nn_sponge_len * Agrif_irhox()196 z1_ispongearea = 1._wp / REAL( ispongearea )197 jspongearea = nn_sponge_len * Agrif_irhoy()198 z1_jspongearea = 1._wp / REAL( jspongearea )199 200 ztabramp(:,:) = 0._wp201 202 ! Trick to remove sponge in 2DV domains:203 IF ( nbcellsx <= 3 ) ispongearea = -1204 IF ( nbcellsy <= 3 ) jspongearea = -1205 206 ! --- West --- !207 IF( lk_west ) THEN208 ind1 = nn_hls + 1 + nbghostcells ! halo + land + nbghostcells209 ind2 = nn_hls + 1 + nbghostcells + ispongearea210 DO ji = mi0(ind1), mi1(ind2)211 DO jj = 1, jpj212 ztabramp(ji,jj) = REAL( ind2 - mig(ji) ) * z1_ispongearea * zmskwest(jj)213 END DO214 END DO215 216 ! ghost cells:217 ind1 = 1218 ind2 = nn_hls + 1 + nbghostcells ! halo + land + nbghostcells219 DO ji = mi0(ind1), mi1(ind2)220 DO jj = 1, jpj221 ztabramp(ji,jj) = zmskwest(jj)222 END DO223 END DO224 ENDIF225 226 ! --- East --- !227 IF( lk_east ) THEN228 ind1 = jpiglo - ( nn_hls + nbghostcells ) - ispongearea229 ind2 = jpiglo - ( nn_hls + nbghostcells ) ! halo + land + nbghostcells - 1230 DO ji = mi0(ind1), mi1(ind2)231 DO jj = 1, jpj232 ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( mig(ji) - ind1 ) * z1_ispongearea) * zmskeast(jj)233 ENDDO234 END DO235 236 ! ghost cells:237 ind1 = jpiglo - ( nn_hls + nbghostcells ) ! halo + land + nbghostcells - 1238 ind2 = jpiglo239 DO ji = mi0(ind1), mi1(ind2)240 DO jj = 1, jpj241 ztabramp(ji,jj) = zmskeast(jj)242 ENDDO243 END DO244 ENDIF245 246 ! --- South --- !247 IF( lk_south ) THEN248 ind1 = nn_hls + 1 + nbghostcells ! halo + land + nbghostcells249 ind2 = nn_hls + 1 + nbghostcells + jspongearea250 DO jj = mj0(ind1), mj1(ind2)251 DO ji = 1, jpi252 ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( ind2 - mjg(jj) ) * z1_jspongearea) * zmsksouth(ji)253 END DO254 END DO255 256 ! ghost cells:257 ind1 = 1258 ind2 = nn_hls + 1 + nbghostcells ! halo + land + nbghostcells259 DO jj = mj0(ind1), mj1(ind2)260 DO ji = 1, jpi261 ztabramp(ji,jj) = zmsksouth(ji)262 END DO263 END DO264 ENDIF265 266 ! --- North --- !267 IF( lk_north ) THEN268 ind1 = jpjglo - ( nn_hls + nbghostcells ) - jspongearea269 ind2 = jpjglo - ( nn_hls + nbghostcells ) ! halo + land + nbghostcells - 1270 DO jj = mj0(ind1), mj1(ind2)271 DO ji = 1, jpi272 ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( mjg(jj) - ind1 ) * z1_jspongearea) * zmsknorth(ji)273 END DO274 END DO275 276 ! ghost cells:277 ind1 = jpjglo - ( nn_hls + nbghostcells ) ! halo + land + nbghostcells - 1278 ind2 = jpjglo279 DO jj = mj0(ind1), mj1(ind2)280 DO ji = 1, jpi281 ztabramp(ji,jj) = zmsknorth(ji)282 END DO283 END DO284 ENDIF285 286 ENDIF287 288 ! Tracers289 IF( .NOT. spongedoneT ) THEN290 fspu(:,:) = 0._wp291 fspv(:,:) = 0._wp292 DO_2D_00_00293 fspu(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji+1,jj ) ) * ssumask(ji,jj)294 fspv(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji ,jj+1) ) * ssvmask(ji,jj)295 END_2D296 ENDIF297 298 ! Dynamics299 IF( .NOT. spongedoneU ) THEN300 fspt(:,:) = 0._wp301 fspf(:,:) = 0._wp302 DO_2D_00_00303 fspt(ji,jj) = ztabramp(ji,jj) * ssmask(ji,jj)304 fspf(ji,jj) = 0.25_wp * ( ztabramp(ji ,jj ) + ztabramp(ji ,jj+1) &305 & +ztabramp(ji+1,jj+1) + ztabramp(ji+1,jj ) ) &306 & * ssvmask(ji,jj) * ssvmask(ji,jj+1)307 END_2D308 ENDIF309 310 IF( .NOT. spongedoneT .AND. .NOT. spongedoneU ) THEN311 CALL lbc_lnk_multi( 'agrif_Sponge', fspu, 'U', 1., fspv, 'V', 1., fspt, 'T', 1., fspf, 'F', 1. )312 spongedoneT = .TRUE.313 spongedoneU = .TRUE.314 ENDIF315 IF( .NOT. spongedoneT ) THEN316 CALL lbc_lnk_multi( 'agrif_Sponge', fspu, 'U', 1., fspv, 'V', 1. )317 spongedoneT = .TRUE.318 ENDIF319 IF( .NOT. spongedoneT .AND. .NOT. spongedoneU ) THEN320 CALL lbc_lnk_multi( 'agrif_Sponge', fspt, 'T', 1., fspf, 'F', 1. )321 spongedoneU = .TRUE.322 ENDIF323 324 #if defined key_vertical325 ! Remove vertical interpolation where not needed:326 DO_2D_00_00327 IF ((fspu(ji-1,jj)==0._wp).AND.(fspu(ji,jj)==0._wp).AND. &328 & (fspv(ji,jj-1)==0._wp).AND.(fspv(ji,jj)==0._wp)) mbkt_parent(ji,jj) = 0329 !330 IF ((fspt(ji+1,jj)==0._wp).AND.(fspt(ji,jj)==0._wp).AND. &331 & (fspf(ji,jj-1)==0._wp).AND.(fspf(ji,jj)==0._wp)) mbku_parent(ji,jj) = 0332 !333 IF ((fspt(ji,jj+1)==0._wp).AND.(fspt(ji,jj)==0._wp).AND. &334 & (fspf(ji-1,jj)==0._wp).AND.(fspf(ji,jj)==0._wp)) mbkv_parent(ji,jj) = 0335 !336 IF ( ssmask(ji,jj) == 0._wp) mbkt_parent(ji,jj) = 0337 IF (ssumask(ji,jj) == 0._wp) mbku_parent(ji,jj) = 0338 IF (ssvmask(ji,jj) == 0._wp) mbkv_parent(ji,jj) = 0339 END_2D340 !341 ztabramp (:,:) = REAL( mbkt_parent (:,:), wp )342 ztabrampu(:,:) = REAL( mbku_parentu(:,:), wp )343 ztabrampv(:,:) = REAL( mbkv_parentv(:,:), wp )344 CALL lbc_lnk_multi( 'Agrif_Sponge', ztabramp, 'T', 1., ztabrampu, 'U', 1., ztabrampv, 'V', 1. )345 mbkt_parent(:,:) = NINT( ztabramp (:,:) )346 mbku_parent(:,:) = NINT( ztabrampu(:,:) )347 mbkv_parent(:,:) = NINT( ztabrampv(:,:) )348 #endif349 !350 #endif351 !352 END SUBROUTINE Agrif_Sponge353 354 SUBROUTINE interptsn_sponge( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )355 !!----------------------------------------------------------------------356 !! *** ROUTINE interptsn_sponge ***357 !!----------------------------------------------------------------------358 INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2359 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres360 LOGICAL , INTENT(in ) :: before361 !362 INTEGER :: ji, jj, jk, jn ! dummy loop indices363 INTEGER :: iku, ikv364 REAL(wp) :: ztsa, zabe1, zabe2, zbtr, zhtot, ztrelax365 REAL(wp), DIMENSION(i1:i2,j1:j2,jpk) :: ztu, ztv366 REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tsbdiff367 ! vertical interpolation:368 REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tabres_child369 REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin370 REAL(wp), DIMENSION(k1:k2) :: h_in371 REAL(wp), DIMENSION(1:jpk) :: h_out372 INTEGER :: N_in, N_out373 !!----------------------------------------------------------------------374 !375 IF( before ) THEN376 DO jn = 1, jpts377 DO jk=k1,k2378 DO jj=j1,j2379 DO ji=i1,i2380 tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kbb_a)381 END DO382 END DO383 END DO384 END DO385 386 # if defined key_vertical387 ! Interpolate thicknesses388 ! Warning: these are masked, hence extrapolated prior interpolation.389 DO jk=k1,k2390 DO jj=j1,j2391 DO ji=i1,i2392 tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kbb_a)393 END DO394 END DO395 END DO396 397 ! Extrapolate thicknesses in partial bottom cells:398 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on399 IF (ln_zps) THEN400 DO jj=j1,j2401 DO ji=i1,i2402 jk = mbkt(ji,jj)403 tabres(ji,jj,jk,jpts+1) = 0._wp404 END DO405 END DO406 END IF407 408 ! Save ssh at last level:409 IF (.NOT.ln_linssh) THEN410 tabres(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kbb_a)*tmask(i1:i2,j1:j2,1)411 ELSE412 tabres(i1:i2,j1:j2,k2,jpts+1) = 0._wp413 END IF414 # endif415 416 ELSE417 !418 # if defined key_vertical419 420 IF (ln_linssh) tabres(i1:i2,j1:j2,k2,n2) = 0._wp421 422 DO jj=j1,j2423 DO ji=i1,i2424 tabres_child(ji,jj,:,:) = 0._wp425 N_in = mbkt_parent(ji,jj)426 zhtot = 0._wp427 DO jk=1,N_in !k2 = jpk of parent grid428 IF (jk==N_in) THEN429 h_in(jk) = ht0_parent(ji,jj) + tabres(ji,jj,k2,n2) - zhtot430 ELSE431 h_in(jk) = tabres(ji,jj,jk,n2)432 ENDIF433 zhtot = zhtot + h_in(jk)434 tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)435 END DO436 N_out = 0437 DO jk=1,jpk ! jpk of child grid438 IF (tmask(ji,jj,jk) == 0) EXIT439 N_out = N_out + 1440 h_out(jk) = e3t(ji,jj,jk,Kbb_a) !Child grid scale factors. Could multiply by e1e2t here instead of division above441 ENDDO442 443 ! Account for small differences in free-surface444 IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN445 h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) )446 ELSE447 h_in(1) = h_in(1) - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) )448 ENDIF449 IF (N_in*N_out > 0) THEN450 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)451 ENDIF452 ENDDO453 ENDDO454 # endif455 456 DO jj=j1,j2457 DO ji=i1,i2458 DO jk=1,jpkm1459 # if defined key_vertical460 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)461 # else462 tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres(ji,jj,jk,1:jpts))*tmask(ji,jj,jk)463 # endif464 ENDDO465 ENDDO466 ENDDO467 468 !* set relaxation time scale469 IF( l_1st_euler .AND. lk_agrif_fstep ) THEN ; ztrelax = rn_trelax_tra / ( rn_Dt )470 ELSE ; ztrelax = rn_trelax_tra / (2._wp * rn_Dt )471 ENDIF472 473 DO jn = 1, jpts474 DO jk = 1, jpkm1475 ztu(i1:i2,j1:j2,jk) = 0._wp476 DO jj = j1,j2477 DO ji = i1,i2-1478 zabe1 = rn_sponge_tra * fspu(ji,jj) * umask(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm_a)479 ztu(ji,jj,jk) = zabe1 * ( tsbdiff(ji+1,jj ,jk,jn) - tsbdiff(ji,jj,jk,jn) )480 END DO481 END DO482 ztv(i1:i2,j1:j2,jk) = 0._wp483 DO ji = i1,i2484 DO jj = j1,j2-1485 zabe2 = rn_sponge_tra * fspv(ji,jj) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a)486 ztv(ji,jj,jk) = zabe2 * ( tsbdiff(ji ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) )487 END DO488 END DO489 !490 IF( ln_zps ) THEN ! set gradient at partial step level491 DO jj = j1,j2492 DO ji = i1,i2493 ! last level494 iku = mbku(ji,jj)495 ikv = mbkv(ji,jj)496 IF( iku == jk ) ztu(ji,jj,jk) = 0._wp497 IF( ikv == jk ) ztv(ji,jj,jk) = 0._wp498 END DO499 END DO500 ENDIF501 END DO502 !503 DO jk = 1, jpkm1504 DO jj = j1+1,j2-1505 DO ji = i1+1,i2-1506 IF (.NOT. tabspongedone_tsn(ji,jj)) THEN507 zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm_a)508 ! horizontal diffusive trends509 ztsa = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) &510 & - ztrelax * fspt(ji,jj) * tsbdiff(ji,jj,jk,jn)511 ! add it to the general tracer trends512 ts(ji,jj,jk,jn,Krhs_a) = ts(ji,jj,jk,jn,Krhs_a) + ztsa513 ENDIF514 END DO515 END DO516 END DO517 !518 END DO519 !520 tabspongedone_tsn(i1+1:i2-1,j1+1:j2-1) = .TRUE.521 !522 ENDIF523 !524 END SUBROUTINE interptsn_sponge525 526 SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before)527 !!---------------------------------------------528 !! *** ROUTINE interpun_sponge ***529 !!---------------------------------------------530 INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2531 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres532 LOGICAL, INTENT(in) :: before533 534 INTEGER :: ji,jj,jk,jmax535 INTEGER :: ind1536 ! sponge parameters537 REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, zhtot, ztrelax538 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ubdiff539 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff540 ! vertical interpolation:541 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child542 REAL(wp), DIMENSION(k1:k2) :: tabin, h_in543 REAL(wp), DIMENSION(1:jpk) :: h_out544 INTEGER ::N_in, N_out545 !!---------------------------------------------546 !547 IF( before ) THEN548 DO jk=k1,k2549 DO jj=j1,j2550 DO ji=i1,i2551 tabres(ji,jj,jk,m1) = uu(ji,jj,jk,Kbb_a)552 # if defined key_vertical553 tabres(ji,jj,jk,m2) = e3u(ji,jj,jk,Kbb_a)*umask(ji,jj,jk)554 # endif555 END DO556 END DO557 END DO558 559 # if defined key_vertical560 ! Extrapolate thicknesses in partial bottom cells:561 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on562 IF (ln_zps) THEN563 DO jj=j1,j2564 DO ji=i1,i2565 jk = mbku(ji,jj)566 tabres(ji,jj,jk,m2) = 0._wp567 END DO568 END DO569 END IF570 ! Save ssh at last level:571 tabres(i1:i2,j1:j2,k2,m2) = 0._wp572 IF (.NOT.ln_linssh) THEN573 ! This vertical sum below should be replaced by the sea-level at U-points (optimization):574 DO jk=1,jpk575 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)576 END DO577 tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hu_0(i1:i2,j1:j2)578 END IF579 # endif580 581 ELSE582 583 # if defined key_vertical584 IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp585 586 DO jj=j1,j2587 DO ji=i1,i2588 tabres_child(ji,jj,:) = 0._wp589 N_in = mbku_parent(ji,jj)590 zhtot = 0._wp591 DO jk=1,N_in592 IF (jk==N_in) THEN593 h_in(jk) = hu0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot594 ELSE595 h_in(jk) = tabres(ji,jj,jk,m2)596 ENDIF597 zhtot = zhtot + h_in(jk)598 tabin(jk) = tabres(ji,jj,jk,m1)599 ENDDO600 !601 N_out = 0602 DO jk=1,jpk603 IF (umask(ji,jj,jk) == 0) EXIT604 N_out = N_out + 1605 h_out(N_out) = e3u(ji,jj,jk,Kbb_a)606 ENDDO607 608 ! Account for small differences in free-surface609 IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN610 h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) )611 ELSE612 h_in(1) = h_in(1) - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) )613 ENDIF614 615 IF (N_in * N_out > 0) THEN616 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)617 ENDIF618 ENDDO619 ENDDO620 621 ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:)622 #else623 ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:)624 #endif625 !* set relaxation time scale626 IF( l_1st_euler .AND. lk_agrif_fstep ) THEN ; ztrelax = rn_trelax_dyn / ( rn_Dt )627 ELSE ; ztrelax = rn_trelax_dyn / (2._wp * rn_Dt )628 ENDIF629 !630 DO jk = 1, jpkm1 ! Horizontal slab631 ! ! ===============632 633 ! ! --------634 ! Horizontal divergence ! div635 ! ! --------636 DO jj = j1,j2637 DO ji = i1+1,i2 ! vector opt.638 zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb_a) * rn_sponge_dyn * fspt(ji,jj)639 hdivdiff(ji,jj,jk) = ( e2u(ji ,jj)*e3u(ji ,jj,jk,Kbb_a) * ubdiff(ji ,jj,jk) &640 & -e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kbb_a) * ubdiff(ji-1,jj,jk) ) * zbtr641 END DO642 END DO643 644 DO jj = j1,j2-1645 DO ji = i1,i2 ! vector opt.646 zbtr = r1_e1e2f(ji,jj) * e3f(ji,jj,jk) * rn_sponge_dyn * fspf(ji,jj)647 rotdiff(ji,jj,jk) = ( -e1u(ji,jj+1) * ubdiff(ji,jj+1,jk) &648 & +e1u(ji,jj ) * ubdiff(ji,jj ,jk) ) * fmask(ji,jj,jk) * zbtr649 END DO650 END DO651 END DO652 !653 DO jj = j1+1, j2-1654 DO ji = i1+1, i2-1 ! vector opt.655 656 IF (.NOT. tabspongedone_u(ji,jj)) THEN657 DO jk = 1, jpkm1 ! Horizontal slab658 ze2u = rotdiff (ji,jj,jk)659 ze1v = hdivdiff(ji,jj,jk)660 ! horizontal diffusive trends661 zua = - ( ze2u - rotdiff (ji,jj-1,jk) ) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) ) &662 & + ( hdivdiff(ji+1,jj,jk) - ze1v ) * r1_e1u(ji,jj) &663 & - ztrelax * fspu(ji,jj) * ubdiff(ji,jj,jk)664 665 ! add it to the general momentum trends666 uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a) + zua667 END DO668 ENDIF669 670 END DO671 END DO672 673 tabspongedone_u(i1+1:i2-1,j1+1:j2-1) = .TRUE.674 675 jmax = j2-1676 ind1 = jpjglo - ( nn_hls + nbghostcells + 2 ) ! North677 DO jj = mj0(ind1), mj1(ind1)678 jmax = MIN(jmax,jj)679 END DO680 681 DO jj = j1+1, jmax682 DO ji = i1+1, i2 ! vector opt.683 684 IF (.NOT. tabspongedone_v(ji,jj)) THEN685 DO jk = 1, jpkm1 ! Horizontal slab686 ze2u = rotdiff (ji,jj,jk)687 ze1v = hdivdiff(ji,jj,jk)688 689 ! horizontal diffusive trends690 zva = + ( ze2u - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) ) &691 + ( hdivdiff(ji,jj+1,jk) - ze1v ) * r1_e2v(ji,jj)692 693 ! add it to the general momentum trends694 vv(ji,jj,jk,Krhs_a) = vv(ji,jj,jk,Krhs_a) + zva695 END DO696 ENDIF697 !698 END DO699 END DO700 !701 tabspongedone_v(i1+1:i2,j1+1:jmax) = .TRUE.702 !703 ENDIF704 !705 END SUBROUTINE interpun_sponge706 707 SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before,nb,ndir)708 !!---------------------------------------------709 !! *** ROUTINE interpvn_sponge ***710 !!---------------------------------------------711 INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2712 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres713 LOGICAL, INTENT(in) :: before714 INTEGER, INTENT(in) :: nb , ndir715 !716 INTEGER :: ji, jj, jk, imax717 INTEGER :: ind1718 ! sponge parameters719 REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, zhtot, ztrelax720 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: vbdiff721 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff722 ! vertical interpolation:723 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child724 REAL(wp), DIMENSION(k1:k2) :: tabin, h_in725 REAL(wp), DIMENSION(1:jpk) :: h_out726 INTEGER :: N_in, N_out727 !!---------------------------------------------728 729 IF( before ) THEN730 DO jk=k1,k2731 DO jj=j1,j2732 34 DO ji=i1,i2 733 35 tabres(ji,jj,jk,m1) = vv(ji,jj,jk,Kbb_a) … … 778 80 zhtot = zhtot + h_in(jk) 779 81 tabin(jk) = tabres(ji,jj,jk,m1) 780 END DO82 END DO 781 83 ! 782 84 N_out = 0 … … 785 87 N_out = N_out + 1 786 88 h_out(N_out) = e3v(ji,jj,jk,Kbb_a) 787 END DO89 END DO 788 90 789 91 ! Account for small differences in free-surface … … 797 99 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) 798 100 ENDIF 799 END DO800 END DO101 END DO 102 END DO 801 103 802 104 vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:) … … 804 106 vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:) 805 107 # endif 806 !* set relaxation time scale807 IF( l_1st_euler .AND. lk_agrif_fstep ) THEN ; ztrelax = rn_trelax_dyn / ( rn_Dt )808 ELSE ; ztrelax = rn_trelax_dyn / (2._wp * rn_Dt )809 ENDIF810 108 ! 811 109 DO jk = 1, jpkm1 ! Horizontal slab … … 817 115 DO jj = j1+1,j2 818 116 DO ji = i1,i2 ! vector opt. 819 zbtr = r 1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb_a) * rn_sponge_dyn * fspt(ji,jj)117 zbtr = rn_sponge_dyn * r1_Dt * fspt(ji,jj) / e3t(ji,jj,jk,Kbb_a) 820 118 hdivdiff(ji,jj,jk) = ( e1v(ji,jj ) * e3v(ji,jj ,jk,Kbb_a) * vbdiff(ji,jj ,jk) & 821 119 & -e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kbb_a) * vbdiff(ji,jj-1,jk) ) * zbtr … … 824 122 DO jj = j1,j2 825 123 DO ji = i1,i2-1 ! vector opt. 826 zbtr = r 1_e1e2f(ji,jj) * e3f(ji,jj,jk) * rn_sponge_dyn * fspf(ji,jj)124 zbtr = rn_sponge_dyn * r1_Dt * fspf(ji,jj) * e3f(ji,jj,jk) 827 125 rotdiff(ji,jj,jk) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj,jk) & 828 126 & -e2v(ji ,jj) * vbdiff(ji ,jj,jk) ) * fmask(ji,jj,jk) * zbtr … … 844 142 IF( .NOT. tabspongedone_u(ji,jj) ) THEN 845 143 DO jk = 1, jpkm1 846 uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a) 144 uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a) & 847 145 & - ( rotdiff (ji ,jj,jk) - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) ) & 848 146 & + ( hdivdiff(ji+1,jj,jk) - hdivdiff(ji,jj ,jk)) * r1_e1u(ji,jj) … … 858 156 IF( .NOT. tabspongedone_v(ji,jj) ) THEN 859 157 DO jk = 1, jpkm1 860 vv(ji,jj,jk,Krhs_a) = vv(ji,jj,jk,Krhs_a) 158 vv(ji,jj,jk,Krhs_a) = vv(ji,jj,jk,Krhs_a) & 861 159 & + ( rotdiff (ji,jj ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) ) & 862 & + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji ,jj,jk) ) * r1_e2v(ji,jj) &863 & - ztrelax* fspv(ji,jj) * vbdiff(ji,jj,jk)160 & + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji ,jj,jk) ) * r1_e2v(ji,jj) & 161 & - rn_trelax_dyn * r1_Dt * fspv(ji,jj) * vbdiff(ji,jj,jk) 864 162 END DO 865 163 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.