- Timestamp:
- 2019-12-09T13:55:34+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_AGRIF-01-05_merged/src/NST/agrif_oce_sponge.F90
r10425 r12123 22 22 USE agrif_oce 23 23 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 24 USE iom 25 USE vremap 24 26 25 27 IMPLICIT NONE … … 58 60 #endif 59 61 ! 62 CALL iom_put( 'agrif_spu', fspu(:,:)) 63 CALL iom_put( 'agrif_spv', fspv(:,:)) 64 ! 60 65 END SUBROUTINE Agrif_Sponge_Tra 61 66 … … 85 90 #endif 86 91 ! 92 CALL iom_put( 'agrif_spt', fspt(:,:)) 93 CALL iom_put( 'agrif_spf', fspf(:,:)) 94 ! 87 95 END SUBROUTINE Agrif_Sponge_dyn 88 96 … … 93 101 !!---------------------------------------------------------------------- 94 102 INTEGER :: ji, jj, ind1, ind2 95 INTEGER :: ispongearea 96 REAL(wp) :: z1_ spongearea103 INTEGER :: ispongearea, jspongearea 104 REAL(wp) :: z1_ispongearea, z1_jspongearea 97 105 REAL(wp), DIMENSION(jpi,jpj) :: ztabramp 98 !!---------------------------------------------------------------------- 99 ! 106 REAL(wp), DIMENSION(jpjmax) :: zmskwest, zmskeast 107 REAL(wp), DIMENSION(jpimax) :: zmsknorth, zmsksouth 108 !!---------------------------------------------------------------------- 109 ! 110 ! Sponge 1d example with: 111 ! iraf = 3 ; nbghost = 3 ; nn_sponge_len = 2 112 ! 113 !coarse : U T U T U T U 114 !| | | | | 115 !fine : t u t u t u t u t u t u t u t u t u t u t 116 !sponge val:0 0 0 1 5/6 4/6 3/6 2/6 1/6 0 0 117 ! | ghost | <-- sponge area -- > | 118 ! | points | | 119 ! |--> dynamical interface 120 100 121 #if defined SPONGE || defined SPONGE_TOP 101 122 IF (( .NOT. spongedoneT ).OR.( .NOT. spongedoneU )) THEN 123 ! 124 ! Retrieve masks at open boundaries: 125 126 ! --- West --- ! 127 ztabramp(:,:) = 0._wp 128 ind1 = 1+nbghostcells 129 DO ji = mi0(ind1), mi1(ind1) 130 ztabramp(ji,:) = ssumask(ji,:) 131 END DO 132 ! 133 zmskwest(:) = 0._wp 134 zmskwest(1:jpj) = MAXVAL(ztabramp(:,:), dim=1) 135 136 ! --- East --- ! 137 ztabramp(:,:) = 0._wp 138 ind1 = jpiglo - nbghostcells - 1 139 DO ji = mi0(ind1), mi1(ind1) 140 ztabramp(ji,:) = ssumask(ji,:) 141 END DO 142 ! 143 zmskeast(:) = 0._wp 144 zmskeast(1:jpj) = MAXVAL(ztabramp(:,:), dim=1) 145 146 ! --- South --- ! 147 ztabramp(:,:) = 0._wp 148 ind1 = 1+nbghostcells 149 DO jj = mj0(ind1), mj1(ind1) 150 ztabramp(:,jj) = ssvmask(:,jj) 151 END DO 152 ! 153 zmsksouth(:) = 0._wp 154 zmsksouth(1:jpi) = MAXVAL(ztabramp(:,:), dim=2) 155 156 ! --- North --- ! 157 ztabramp(:,:) = 0._wp 158 ind1 = jpjglo - nbghostcells - 1 159 DO jj = mj0(ind1), mj1(ind1) 160 ztabramp(:,jj) = ssvmask(:,jj) 161 END DO 162 ! 163 zmsknorth(:) = 0._wp 164 zmsknorth(1:jpi) = MAXVAL(ztabramp(:,:), dim=2) 165 166 #if defined key_mpp_mpi 167 CALL mpp_max( 'AGRIF_sponge', zmskwest(:) , jpjmax ) 168 CALL mpp_max( 'AGRIF_Sponge', zmskeast(:) , jpjmax ) 169 CALL mpp_max( 'AGRIF_Sponge', zmsksouth(:), jpimax ) 170 CALL mpp_max( 'AGRIF_Sponge', zmsknorth(:), jpimax ) 171 #endif 172 102 173 ! Define ramp from boundaries towards domain interior at T-points 103 174 ! Store it in ztabramp 104 175 105 ispongearea = 1 + nn_sponge_len * Agrif_irhox() 106 z1_spongearea = 1._wp / REAL( ispongearea ) 176 ispongearea = nn_sponge_len * Agrif_irhox() 177 z1_ispongearea = 1._wp / REAL( ispongearea ) 178 jspongearea = nn_sponge_len * Agrif_irhoy() 179 z1_jspongearea = 1._wp / REAL( jspongearea ) 107 180 108 181 ztabramp(:,:) = 0._wp 109 182 183 ! Trick to remove sponge in 2DV domains: 184 IF ( nbcellsx <= 3 ) ispongearea = -1 185 IF ( nbcellsy <= 3 ) jspongearea = -1 186 110 187 ! --- West --- ! 111 IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN 112 ind1 = 1+nbghostcells 113 ind2 = 1+nbghostcells + ispongearea 188 ind1 = 1+nbghostcells 189 ind2 = 1+nbghostcells + ispongearea 190 DO ji = mi0(ind1), mi1(ind2) 191 DO jj = 1, jpj 192 ztabramp(ji,jj) = REAL( ind2 - mig(ji) ) * z1_ispongearea * zmskwest(jj) 193 END DO 194 END DO 195 196 ! ghost cells: 197 ind1 = 1 198 ind2 = nbghostcells + 1 199 DO ji = mi0(ind1), mi1(ind2) 200 DO jj = 1, jpj 201 ztabramp(ji,jj) = zmskwest(jj) 202 END DO 203 END DO 204 205 ! --- East --- ! 206 ind1 = jpiglo - nbghostcells - ispongearea 207 ind2 = jpiglo - nbghostcells 208 DO ji = mi0(ind1), mi1(ind2) 114 209 DO jj = 1, jpj 115 DO ji = ind1, ind2 116 ztabramp(ji,jj) = REAL( ind2 - ji ) * z1_spongearea * umask(ind1,jj,1) 117 END DO 210 ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( mig(ji) - ind1 ) * z1_ispongearea) * zmskeast(jj) 118 211 ENDDO 119 END IF120 121 ! --- East --- !122 IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN123 ind1 = nlci - nbghostcells - ispongearea124 ind2 = nlci - nbghostcells212 END DO 213 214 ! ghost cells: 215 ind1 = jpiglo - nbghostcells 216 ind2 = jpiglo 217 DO ji = mi0(ind1), mi1(ind2) 125 218 DO jj = 1, jpj 126 DO ji = ind1, ind2 127 ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( ji - ind1 ) * z1_spongearea * umask(ind2-1,jj,1) ) 128 ENDDO 219 ztabramp(ji,jj) = zmskeast(jj) 129 220 ENDDO 130 END IF221 END DO 131 222 132 223 ! --- South --- ! 133 IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN 134 ind1 = 1+nbghostcells 135 ind2 = 1+nbghostcells + ispongearea 136 DO jj = ind1, ind2 137 DO ji = 1, jpi 138 ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( ind2 - jj ) * z1_spongearea * vmask(ji,ind1,1) ) 139 END DO 140 ENDDO 141 ENDIF 224 ind1 = 1+nbghostcells 225 ind2 = 1+nbghostcells + jspongearea 226 DO jj = mj0(ind1), mj1(ind2) 227 DO ji = 1, jpi 228 ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( ind2 - mjg(jj) ) * z1_jspongearea) * zmsksouth(ji) 229 END DO 230 END DO 231 232 ! ghost cells: 233 ind1 = 1 234 ind2 = nbghostcells + 1 235 DO jj = mj0(ind1), mj1(ind2) 236 DO ji = 1, jpi 237 ztabramp(ji,jj) = zmsksouth(ji) 238 END DO 239 END DO 142 240 143 241 ! --- North --- ! 144 IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN 145 ind1 = nlcj - nbghostcells - ispongearea 146 ind2 = nlcj - nbghostcells 147 DO jj = ind1, ind2 148 DO ji = 1, jpi 149 ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( jj - ind1 ) * z1_spongearea * vmask(ji,ind2-1,1) ) 150 END DO 151 ENDDO 152 ENDIF 242 ind1 = jpjglo - nbghostcells - jspongearea 243 ind2 = jpjglo - nbghostcells 244 DO jj = mj0(ind1), mj1(ind2) 245 DO ji = 1, jpi 246 ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( mjg(jj) - ind1 ) * z1_jspongearea) * zmsknorth(ji) 247 END DO 248 END DO 249 250 ! ghost cells: 251 ind1 = jpjglo - nbghostcells 252 ind2 = jpjglo 253 DO jj = mj0(ind1), mj1(ind2) 254 DO ji = 1, jpi 255 ztabramp(ji,jj) = zmsknorth(ji) 256 END DO 257 END DO 153 258 154 259 ENDIF … … 156 261 ! Tracers 157 262 IF( .NOT. spongedoneT ) THEN 158 fs aht_spu(:,:) = 0._wp159 fs aht_spv(:,:) = 0._wp263 fspu(:,:) = 0._wp 264 fspv(:,:) = 0._wp 160 265 DO jj = 2, jpjm1 161 266 DO ji = 2, jpim1 ! vector opt. 162 fs aht_spu(ji,jj) = 0.5_wp * visc_tra * ( ztabramp(ji,jj) + ztabramp(ji+1,jj ))163 fs aht_spv(ji,jj) = 0.5_wp * visc_tra * ( ztabramp(ji,jj) + ztabramp(ji ,jj+1))164 END DO 165 END DO 166 CALL lbc_lnk( 'agrif_ oce_sponge', fsaht_spu, 'U', 1. ) ! Lateral boundary conditions167 CALL lbc_lnk( 'agrif_ oce_sponge', fsaht_spv, 'V', 1. )168 267 fspu(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji+1,jj ) ) * ssumask(ji,jj) 268 fspv(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji ,jj+1) ) * ssvmask(ji,jj) 269 END DO 270 END DO 271 CALL lbc_lnk( 'agrif_Sponge', fspu, 'U', 1. ) ! Lateral boundary conditions 272 CALL lbc_lnk( 'agrif_Sponge', fspv, 'V', 1. ) 273 169 274 spongedoneT = .TRUE. 170 275 ENDIF … … 172 277 ! Dynamics 173 278 IF( .NOT. spongedoneU ) THEN 174 fs ahm_spt(:,:) = 0._wp175 fs ahm_spf(:,:) = 0._wp279 fspt(:,:) = 0._wp 280 fspf(:,:) = 0._wp 176 281 DO jj = 2, jpjm1 177 282 DO ji = 2, jpim1 ! vector opt. 178 fsahm_spt(ji,jj) = visc_dyn * ztabramp(ji,jj) 179 fsahm_spf(ji,jj) = 0.25_wp * visc_dyn * ( ztabramp(ji ,jj ) + ztabramp(ji ,jj+1) & 180 & +ztabramp(ji+1,jj+1) + ztabramp(ji+1,jj ) ) 181 END DO 182 END DO 183 CALL lbc_lnk( 'agrif_oce_sponge', fsahm_spt, 'T', 1. ) ! Lateral boundary conditions 184 CALL lbc_lnk( 'agrif_oce_sponge', fsahm_spf, 'F', 1. ) 283 fspt(ji,jj) = ztabramp(ji,jj) * ssmask(ji,jj) 284 fspf(ji,jj) = 0.25_wp * ( ztabramp(ji ,jj ) + ztabramp(ji ,jj+1) & 285 & +ztabramp(ji+1,jj+1) + ztabramp(ji+1,jj ) ) & 286 & * ssvmask(ji,jj) * ssvmask(ji,jj+1) 287 END DO 288 END DO 289 CALL lbc_lnk( 'agrif_Sponge', fspt, 'T', 1. ) ! Lateral boundary conditions 290 CALL lbc_lnk( 'agrif_Sponge', fspf, 'F', 1. ) 185 291 186 292 spongedoneU = .TRUE. 187 293 ENDIF 294 295 #if defined key_vertical 296 ! Remove vertical interpolation where not needed: 297 DO jj = 2, jpjm1 298 DO ji = 2, jpim1 299 IF ((fspu(ji-1,jj)==0._wp).AND.(fspu(ji,jj)==0._wp).AND. & 300 & (fspv(ji,jj-1)==0._wp).AND.(fspv(ji,jj)==0._wp)) mbkt_parent(ji,jj) = 0 301 ! 302 IF ((fspt(ji+1,jj)==0._wp).AND.(fspt(ji,jj)==0._wp).AND. & 303 & (fspf(ji,jj-1)==0._wp).AND.(fspf(ji,jj)==0._wp)) mbku_parent(ji,jj) = 0 304 ! 305 IF ((fspt(ji,jj+1)==0._wp).AND.(fspt(ji,jj)==0._wp).AND. & 306 & (fspf(ji-1,jj)==0._wp).AND.(fspf(ji,jj)==0._wp)) mbkv_parent(ji,jj) = 0 307 ! 308 IF ( ssmask(ji,jj) == 0._wp) mbkt_parent(ji,jj) = 0 309 IF (ssumask(ji,jj) == 0._wp) mbku_parent(ji,jj) = 0 310 IF (ssvmask(ji,jj) == 0._wp) mbkv_parent(ji,jj) = 0 311 END DO 312 END DO 313 ! 314 ztabramp(:,:) = REAL( mbkt_parent(:,:), wp ) ; CALL lbc_lnk( 'Agrif_Sponge', ztabramp, 'T', 1. ) 315 mbkt_parent(:,:) = NINT( ztabramp(:,:) ) 316 ztabramp(:,:) = REAL( mbku_parent(:,:), wp ) ; CALL lbc_lnk( 'Agrif_Sponge', ztabramp, 'U', 1. ) 317 mbku_parent(:,:) = NINT( ztabramp(:,:) ) 318 ztabramp(:,:) = REAL( mbkv_parent(:,:), wp ) ; CALL lbc_lnk( 'Agrif_Sponge', ztabramp, 'V', 1. ) 319 mbkv_parent(:,:) = NINT( ztabramp(:,:) ) 320 #endif 188 321 ! 189 322 #endif … … 201 334 INTEGER :: ji, jj, jk, jn ! dummy loop indices 202 335 INTEGER :: iku, ikv 203 REAL(wp) :: ztsa, zabe1, zabe2, zbtr 336 REAL(wp) :: ztsa, zabe1, zabe2, zbtr, zhtot, ztrelax 204 337 REAL(wp), DIMENSION(i1:i2,j1:j2,jpk) :: ztu, ztv 205 338 REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tsbdiff … … 210 343 REAL(wp), DIMENSION(1:jpk) :: h_out 211 344 INTEGER :: N_in, N_out 212 REAL(wp) :: h_diff213 345 !!---------------------------------------------------------------------- 214 346 ! … … 225 357 226 358 # if defined key_vertical 227 DO jk=k1,k2 228 DO jj=j1,j2 229 DO ji=i1,i2 230 tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk) 231 END DO 232 END DO 233 END DO 359 ! Interpolate thicknesses 360 ! Warning: these are masked, hence extrapolated prior interpolation. 361 DO jk=k1,k2 362 DO jj=j1,j2 363 DO ji=i1,i2 364 tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_b(ji,jj,jk) 365 END DO 366 END DO 367 END DO 368 369 ! Extrapolate thicknesses in partial bottom cells: 370 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 371 IF (ln_zps) THEN 372 DO jj=j1,j2 373 DO ji=i1,i2 374 jk = mbkt(ji,jj) 375 tabres(ji,jj,jk,jpts+1) = 0._wp 376 END DO 377 END DO 378 END IF 379 380 ! Save ssh at last level: 381 IF (.NOT.ln_linssh) THEN 382 tabres(i1:i2,j1:j2,k2,jpts+1) = sshb(i1:i2,j1:j2)*tmask(i1:i2,j1:j2,1) 383 ELSE 384 tabres(i1:i2,j1:j2,k2,jpts+1) = 0._wp 385 END IF 234 386 # endif 235 387 … … 237 389 ! 238 390 # if defined key_vertical 239 tabres_child(:,:,:,:) = 0. 391 392 IF (ln_linssh) tabres(i1:i2,j1:j2,k2,n2) = 0._wp 393 240 394 DO jj=j1,j2 241 395 DO ji=i1,i2 242 N_in = 0 243 DO jk=k1,k2 !k2 = jpk of parent grid 244 IF (tabres(ji,jj,jk,n2) == 0) EXIT 245 N_in = N_in + 1 396 tabres_child(ji,jj,:,:) = 0._wp 397 N_in = mbkt_parent(ji,jj) 398 zhtot = 0._wp 399 DO jk=1,N_in !k2 = jpk of parent grid 400 IF (jk==N_in) THEN 401 h_in(jk) = ht0_parent(ji,jj) + tabres(ji,jj,k2,n2) - zhtot 402 ELSE 403 h_in(jk) = tabres(ji,jj,jk,n2) 404 ENDIF 405 zhtot = zhtot + h_in(jk) 246 406 tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1) 247 h_in(N_in) = tabres(ji,jj,jk,n2)248 407 END DO 249 408 N_out = 0 … … 251 410 IF (tmask(ji,jj,jk) == 0) EXIT 252 411 N_out = N_out + 1 253 h_out(jk) = e3t_ n(ji,jj,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above412 h_out(jk) = e3t_b(ji,jj,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above 254 413 ENDDO 255 IF (N_in > 0) THEN 256 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 257 tabres(ji,jj,k2,:) = tabres(ji,jj,k2-1,:) !what is this line for????? 258 DO jn=1,jpts 259 call reconstructandremap(tabin(1:N_in,jn),h_in,tabres_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 260 ENDDO 414 415 ! Account for small differences in free-surface 416 IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN 417 h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) ) 418 ELSE 419 h_in(1) = h_in(1) - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) ) 420 ENDIF 421 IF (N_in*N_out > 0) THEN 422 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) 261 423 ENDIF 262 424 ENDDO … … 268 430 DO jk=1,jpkm1 269 431 # if defined key_vertical 270 tsbdiff(ji,jj,jk,1:jpts) = tsb(ji,jj,jk,1:jpts) - tabres_child(ji,jj,jk,1:jpts)432 tsbdiff(ji,jj,jk,1:jpts) = (tsb(ji,jj,jk,1:jpts) - tabres_child(ji,jj,jk,1:jpts)) * tmask(ji,jj,jk) 271 433 # else 272 tsbdiff(ji,jj,jk,1:jpts) = tsb(ji,jj,jk,1:jpts) - tabres(ji,jj,jk,1:jpts)434 tsbdiff(ji,jj,jk,1:jpts) = (tsb(ji,jj,jk,1:jpts) - tabres(ji,jj,jk,1:jpts))*tmask(ji,jj,jk) 273 435 # endif 274 436 ENDDO 275 437 ENDDO 276 438 ENDDO 439 440 !* set relaxation time scale 441 IF( neuler == 0 .AND. lk_agrif_fstep ) THEN ; ztrelax = rn_trelax_tra / ( rdt ) 442 ELSE ; ztrelax = rn_trelax_tra / (2._wp * rdt ) 443 ENDIF 277 444 278 445 DO jn = 1, jpts … … 281 448 DO jj = j1,j2 282 449 DO ji = i1,i2-1 283 zabe1 = fsaht_spu(ji,jj) * umask(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk)450 zabe1 = rn_sponge_tra * fspu(ji,jj) * umask(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk) 284 451 ztu(ji,jj,jk) = zabe1 * ( tsbdiff(ji+1,jj ,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 285 452 END DO … … 288 455 DO ji = i1,i2 289 456 DO jj = j1,j2-1 290 zabe2 = fsaht_spv(ji,jj) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk)457 zabe2 = rn_sponge_tra * fspv(ji,jj) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk) 291 458 ztv(ji,jj,jk) = zabe2 * ( tsbdiff(ji ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 292 459 END DO … … 312 479 zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 313 480 ! horizontal diffusive trends 314 ztsa = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) 481 ztsa = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) & 482 & - ztrelax * fspt(ji,jj) * tsbdiff(ji,jj,jk,jn) 315 483 ! add it to the general tracer trends 316 484 tsa(ji,jj,jk,jn) = tsa(ji,jj,jk,jn) + ztsa … … 339 507 340 508 ! sponge parameters 341 REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, h_diff509 REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, zhtot, ztrelax 342 510 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ubdiff 343 511 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff … … 346 514 REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 347 515 REAL(wp), DIMENSION(1:jpk) :: h_out 348 INTEGER ::N_in, N_out516 INTEGER ::N_in, N_out 349 517 !!--------------------------------------------- 350 518 ! 351 519 IF( before ) THEN 352 DO jk= 1,jpkm1520 DO jk=k1,k2 353 521 DO jj=j1,j2 354 522 DO ji=i1,i2 355 523 tabres(ji,jj,jk,m1) = ub(ji,jj,jk) 356 524 # if defined key_vertical 357 tabres(ji,jj,jk,m2) = e3u_ n(ji,jj,jk)*umask(ji,jj,jk)525 tabres(ji,jj,jk,m2) = e3u_b(ji,jj,jk)*umask(ji,jj,jk) 358 526 # endif 359 527 END DO … … 361 529 END DO 362 530 531 # if defined key_vertical 532 ! Extrapolate thicknesses in partial bottom cells: 533 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 534 IF (ln_zps) THEN 535 DO jj=j1,j2 536 DO ji=i1,i2 537 jk = mbku(ji,jj) 538 tabres(ji,jj,jk,m2) = 0._wp 539 END DO 540 END DO 541 END IF 542 ! Save ssh at last level: 543 tabres(i1:i2,j1:j2,k2,m2) = 0._wp 544 IF (.NOT.ln_linssh) THEN 545 ! This vertical sum below should be replaced by the sea-level at U-points (optimization): 546 DO jk=1,jpk 547 tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) + e3u_b(i1:i2,j1:j2,jk) * umask(i1:i2,j1:j2,jk) 548 END DO 549 tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hu_0(i1:i2,j1:j2) 550 END IF 551 # endif 552 363 553 ELSE 364 554 365 555 # if defined key_vertical 366 tabres_child(:,:,:) = 0._wp 556 IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp 557 367 558 DO jj=j1,j2 368 559 DO ji=i1,i2 369 N_in = 0 370 DO jk=k1,k2 371 IF (tabres(ji,jj,jk,m2) == 0) EXIT 372 N_in = N_in + 1 560 tabres_child(ji,jj,:) = 0._wp 561 N_in = mbku_parent(ji,jj) 562 zhtot = 0._wp 563 DO jk=1,N_in 564 IF (jk==N_in) THEN 565 h_in(jk) = hu0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot 566 ELSE 567 h_in(jk) = tabres(ji,jj,jk,m2) 568 ENDIF 569 zhtot = zhtot + h_in(jk) 373 570 tabin(jk) = tabres(ji,jj,jk,m1) 374 h_in(N_in) = tabres(ji,jj,jk,m2) 375 ENDDO 376 ! 377 IF (N_in == 0) THEN 378 tabres_child(ji,jj,:) = 0. 379 CYCLE 380 ENDIF 381 382 N_out = 0 383 DO jk=1,jpk 384 if (umask(ji,jj,jk) == 0) EXIT 385 N_out = N_out + 1 386 h_out(N_out) = e3u_n(ji,jj,jk) 387 ENDDO 388 389 IF (N_out == 0) THEN 390 tabres_child(ji,jj,:) = 0. 391 CYCLE 392 ENDIF 393 394 IF (N_in * N_out > 0) THEN 395 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 396 if (h_diff < -1.e4) then 397 print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in)) 398 endif 399 ENDIF 400 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) 401 571 ENDDO 572 ! 573 N_out = 0 574 DO jk=1,jpk 575 IF (umask(ji,jj,jk) == 0) EXIT 576 N_out = N_out + 1 577 h_out(N_out) = e3u_b(ji,jj,jk) 578 ENDDO 579 580 ! Account for small differences in free-surface 581 IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN 582 h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) ) 583 ELSE 584 h_in(1) = h_in(1) - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) ) 585 ENDIF 586 587 IF (N_in * N_out > 0) THEN 588 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) 589 ENDIF 402 590 ENDDO 403 591 ENDDO … … 407 595 ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:) 408 596 #endif 597 !* set relaxation time scale 598 IF( neuler == 0 .AND. lk_agrif_fstep ) THEN ; ztrelax = rn_trelax_dyn / ( rdt ) 599 ELSE ; ztrelax = rn_trelax_dyn / (2._wp * rdt ) 600 ENDIF 409 601 ! 410 602 DO jk = 1, jpkm1 ! Horizontal slab … … 416 608 DO jj = j1,j2 417 609 DO ji = i1+1,i2 ! vector opt. 418 zbtr = r1_e1e2t(ji,jj) / e3t_ n(ji,jj,jk) * fsahm_spt(ji,jj)419 hdivdiff(ji,jj,jk) = ( e2u(ji ,jj)*e3u_ n(ji ,jj,jk) * ubdiff(ji ,jj,jk) &420 & -e2u(ji-1,jj)*e3u_ n(ji-1,jj,jk) * ubdiff(ji-1,jj,jk) ) * zbtr610 zbtr = r1_e1e2t(ji,jj) / e3t_b(ji,jj,jk) * rn_sponge_dyn * fspt(ji,jj) 611 hdivdiff(ji,jj,jk) = ( e2u(ji ,jj)*e3u_b(ji ,jj,jk) * ubdiff(ji ,jj,jk) & 612 & -e2u(ji-1,jj)*e3u_b(ji-1,jj,jk) * ubdiff(ji-1,jj,jk) ) * zbtr 421 613 END DO 422 614 END DO … … 424 616 DO jj = j1,j2-1 425 617 DO ji = i1,i2 ! vector opt. 426 zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * fsahm_spf(ji,jj)618 zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * rn_sponge_dyn * fspf(ji,jj) 427 619 rotdiff(ji,jj,jk) = ( -e1u(ji,jj+1) * ubdiff(ji,jj+1,jk) & 428 620 & +e1u(ji,jj ) * ubdiff(ji,jj ,jk) ) * fmask(ji,jj,jk) * zbtr … … 440 632 ! horizontal diffusive trends 441 633 zua = - ( ze2u - rotdiff (ji,jj-1,jk) ) / ( e2u(ji,jj) * e3u_n(ji,jj,jk) ) & 442 + ( hdivdiff(ji+1,jj,jk) - ze1v ) * r1_e1u(ji,jj) 634 & + ( hdivdiff(ji+1,jj,jk) - ze1v ) * r1_e1u(ji,jj) & 635 & - ztrelax * fspu(ji,jj) * ubdiff(ji,jj,jk) 443 636 444 637 ! add it to the general momentum trends 445 ua(ji,jj,jk) = ua(ji,jj,jk) + zua 446 638 ua(ji,jj,jk) = ua(ji,jj,jk) + zua 447 639 END DO 448 640 ENDIF … … 492 684 ! 493 685 INTEGER :: ji, jj, jk, imax 494 REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, h_diff686 REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, zhtot, ztrelax 495 687 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: vbdiff 496 688 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff … … 503 695 504 696 IF( before ) THEN 505 DO jk= 1,jpkm1697 DO jk=k1,k2 506 698 DO jj=j1,j2 507 699 DO ji=i1,i2 508 700 tabres(ji,jj,jk,m1) = vb(ji,jj,jk) 509 701 # if defined key_vertical 510 tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v_ n(ji,jj,jk)702 tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v_b(ji,jj,jk) 511 703 # endif 512 704 END DO 513 705 END DO 514 706 END DO 707 708 # if defined key_vertical 709 ! Extrapolate thicknesses in partial bottom cells: 710 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 711 IF (ln_zps) THEN 712 DO jj=j1,j2 713 DO ji=i1,i2 714 jk = mbkv(ji,jj) 715 tabres(ji,jj,jk,m2) = 0._wp 716 END DO 717 END DO 718 END IF 719 ! Save ssh at last level: 720 tabres(i1:i2,j1:j2,k2,m2) = 0._wp 721 IF (.NOT.ln_linssh) THEN 722 ! This vertical sum below should be replaced by the sea-level at V-points (optimization): 723 DO jk=1,jpk 724 tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) + e3v_b(i1:i2,j1:j2,jk) * vmask(i1:i2,j1:j2,jk) 725 END DO 726 tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hv_0(i1:i2,j1:j2) 727 END IF 728 # endif 729 515 730 ELSE 516 731 517 732 # if defined key_vertical 518 tabres_child(:,:,:) = 0._wp733 IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp 519 734 DO jj=j1,j2 520 735 DO ji=i1,i2 521 N_in = 0 522 DO jk=k1,k2 523 IF (tabres(ji,jj,jk,m2) == 0) EXIT 524 N_in = N_in + 1 736 tabres_child(ji,jj,:) = 0._wp 737 N_in = mbkv_parent(ji,jj) 738 zhtot = 0._wp 739 DO jk=1,N_in 740 IF (jk==N_in) THEN 741 h_in(jk) = hv0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot 742 ELSE 743 h_in(jk) = tabres(ji,jj,jk,m2) 744 ENDIF 745 zhtot = zhtot + h_in(jk) 525 746 tabin(jk) = tabres(ji,jj,jk,m1) 526 h_in(N_in) = tabres(ji,jj,jk,m2) 527 ENDDO 747 ENDDO 748 ! 749 N_out = 0 750 DO jk=1,jpk 751 IF (vmask(ji,jj,jk) == 0) EXIT 752 N_out = N_out + 1 753 h_out(N_out) = e3v_b(ji,jj,jk) 754 ENDDO 755 756 ! Account for small differences in free-surface 757 IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN 758 h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) ) 759 ELSE 760 h_in(1) = h_in(1) - ( sum(h_in(1:N_in))-sum(h_out(1:N_out)) ) 761 ENDIF 528 762 529 IF (N_in == 0) THEN 530 tabres_child(ji,jj,:) = 0. 531 CYCLE 532 ENDIF 533 534 N_out = 0 535 DO jk=1,jpk 536 if (vmask(ji,jj,jk) == 0) EXIT 537 N_out = N_out + 1 538 h_out(N_out) = e3v_n(ji,jj,jk) 539 ENDDO 540 541 IF (N_in * N_out > 0) THEN 542 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 543 if (h_diff < -1.e4) then 544 print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in)) 545 endif 546 ENDIF 547 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) 763 IF (N_in * N_out > 0) THEN 764 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) 765 ENDIF 548 766 ENDDO 549 767 ENDDO … … 553 771 vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:) 554 772 # endif 773 !* set relaxation time scale 774 IF( neuler == 0 .AND. lk_agrif_fstep ) THEN ; ztrelax = rn_trelax_dyn / ( rdt ) 775 ELSE ; ztrelax = rn_trelax_dyn / (2._wp * rdt ) 776 ENDIF 555 777 ! 556 778 DO jk = 1, jpkm1 ! Horizontal slab … … 562 784 DO jj = j1+1,j2 563 785 DO ji = i1,i2 ! vector opt. 564 zbtr = r1_e1e2t(ji,jj) / e3t_ n(ji,jj,jk) * fsahm_spt(ji,jj)565 hdivdiff(ji,jj,jk) = ( e1v(ji,jj ) * e3v_ n(ji,jj ,jk) * vbdiff(ji,jj ,jk) &566 & -e1v(ji,jj-1) * e3v_ n(ji,jj-1,jk) * vbdiff(ji,jj-1,jk) ) * zbtr786 zbtr = r1_e1e2t(ji,jj) / e3t_b(ji,jj,jk) * rn_sponge_dyn * fspt(ji,jj) 787 hdivdiff(ji,jj,jk) = ( e1v(ji,jj ) * e3v_b(ji,jj ,jk) * vbdiff(ji,jj ,jk) & 788 & -e1v(ji,jj-1) * e3v_b(ji,jj-1,jk) * vbdiff(ji,jj-1,jk) ) * zbtr 567 789 END DO 568 790 END DO 569 791 DO jj = j1,j2 570 792 DO ji = i1,i2-1 ! vector opt. 571 zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * fsahm_spf(ji,jj)793 zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * rn_sponge_dyn * fspf(ji,jj) 572 794 rotdiff(ji,jj,jk) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj,jk) & 573 795 & -e2v(ji ,jj) * vbdiff(ji ,jj,jk) ) * fmask(ji,jj,jk) * zbtr … … 602 824 va(ji,jj,jk) = va(ji,jj,jk) & 603 825 & + ( rotdiff (ji,jj ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v_n(ji,jj,jk) ) & 604 & + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji ,jj,jk) ) * r1_e2v(ji,jj) 826 & + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji ,jj,jk) ) * r1_e2v(ji,jj) & 827 & - ztrelax * fspv(ji,jj) * vbdiff(ji,jj,jk) 605 828 END DO 606 829 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
