- Timestamp:
- 2019-12-11T16:56:06+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/NST/agrif_oce_sponge.F90
r10425 r12191 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 ! JC: SPONGE MASKING TO BE SORTED OUT: 166 zmskwest(:) = 1._wp 167 zmskeast(:) = 1._wp 168 zmsknorth(:) = 1._wp 169 zmsksouth(:) = 1._wp 170 #if defined key_mpp_mpi 171 ! CALL mpp_max( 'AGRIF_sponge', zmskwest(:) , jpjmax ) 172 ! CALL mpp_max( 'AGRIF_Sponge', zmskeast(:) , jpjmax ) 173 ! CALL mpp_max( 'AGRIF_Sponge', zmsksouth(:), jpimax ) 174 ! CALL mpp_max( 'AGRIF_Sponge', zmsknorth(:), jpimax ) 175 #endif 176 102 177 ! Define ramp from boundaries towards domain interior at T-points 103 178 ! Store it in ztabramp 104 179 105 ispongearea = 1 + nn_sponge_len * Agrif_irhox() 106 z1_spongearea = 1._wp / REAL( ispongearea ) 180 ispongearea = nn_sponge_len * Agrif_irhox() 181 z1_ispongearea = 1._wp / REAL( ispongearea ) 182 jspongearea = nn_sponge_len * Agrif_irhoy() 183 z1_jspongearea = 1._wp / REAL( jspongearea ) 107 184 108 185 ztabramp(:,:) = 0._wp 109 186 187 ! Trick to remove sponge in 2DV domains: 188 IF ( nbcellsx <= 3 ) ispongearea = -1 189 IF ( nbcellsy <= 3 ) jspongearea = -1 190 110 191 ! --- West --- ! 111 IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN 112 ind1 = 1+nbghostcells 113 ind2 = 1+nbghostcells + ispongearea 192 ind1 = 1+nbghostcells 193 ind2 = 1+nbghostcells + ispongearea 194 DO ji = mi0(ind1), mi1(ind2) 195 DO jj = 1, jpj 196 ztabramp(ji,jj) = REAL( ind2 - mig(ji) ) * z1_ispongearea * zmskwest(jj) 197 END DO 198 END DO 199 200 ! ghost cells: 201 ind1 = 1 202 ind2 = nbghostcells + 1 203 DO ji = mi0(ind1), mi1(ind2) 204 DO jj = 1, jpj 205 ztabramp(ji,jj) = zmskwest(jj) 206 END DO 207 END DO 208 209 ! --- East --- ! 210 ind1 = jpiglo - nbghostcells - ispongearea 211 ind2 = jpiglo - nbghostcells 212 DO ji = mi0(ind1), mi1(ind2) 114 213 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 214 ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( mig(ji) - ind1 ) * z1_ispongearea) * zmskeast(jj) 118 215 ENDDO 119 END IF120 121 ! --- East --- !122 IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN123 ind1 = nlci - nbghostcells - ispongearea124 ind2 = nlci - nbghostcells216 END DO 217 218 ! ghost cells: 219 ind1 = jpiglo - nbghostcells 220 ind2 = jpiglo 221 DO ji = mi0(ind1), mi1(ind2) 125 222 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 223 ztabramp(ji,jj) = zmskeast(jj) 129 224 ENDDO 130 END IF225 END DO 131 226 132 227 ! --- 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 228 ind1 = 1+nbghostcells 229 ind2 = 1+nbghostcells + jspongearea 230 DO jj = mj0(ind1), mj1(ind2) 231 DO ji = 1, jpi 232 ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( ind2 - mjg(jj) ) * z1_jspongearea) * zmsksouth(ji) 233 END DO 234 END DO 235 236 ! ghost cells: 237 ind1 = 1 238 ind2 = nbghostcells + 1 239 DO jj = mj0(ind1), mj1(ind2) 240 DO ji = 1, jpi 241 ztabramp(ji,jj) = zmsksouth(ji) 242 END DO 243 END DO 142 244 143 245 ! --- 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 246 ind1 = jpjglo - nbghostcells - jspongearea 247 ind2 = jpjglo - nbghostcells 248 DO jj = mj0(ind1), mj1(ind2) 249 DO ji = 1, jpi 250 ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( mjg(jj) - ind1 ) * z1_jspongearea) * zmsknorth(ji) 251 END DO 252 END DO 253 254 ! ghost cells: 255 ind1 = jpjglo - nbghostcells 256 ind2 = jpjglo 257 DO jj = mj0(ind1), mj1(ind2) 258 DO ji = 1, jpi 259 ztabramp(ji,jj) = zmsknorth(ji) 260 END DO 261 END DO 153 262 154 263 ENDIF … … 156 265 ! Tracers 157 266 IF( .NOT. spongedoneT ) THEN 158 fs aht_spu(:,:) = 0._wp159 fs aht_spv(:,:) = 0._wp267 fspu(:,:) = 0._wp 268 fspv(:,:) = 0._wp 160 269 DO jj = 2, jpjm1 161 270 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 271 fspu(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji+1,jj ) ) * ssumask(ji,jj) 272 fspv(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji ,jj+1) ) * ssvmask(ji,jj) 273 END DO 274 END DO 275 CALL lbc_lnk( 'agrif_Sponge', fspu, 'U', 1. ) ! Lateral boundary conditions 276 CALL lbc_lnk( 'agrif_Sponge', fspv, 'V', 1. ) 277 169 278 spongedoneT = .TRUE. 170 279 ENDIF … … 172 281 ! Dynamics 173 282 IF( .NOT. spongedoneU ) THEN 174 fs ahm_spt(:,:) = 0._wp175 fs ahm_spf(:,:) = 0._wp283 fspt(:,:) = 0._wp 284 fspf(:,:) = 0._wp 176 285 DO jj = 2, jpjm1 177 286 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. ) 287 fspt(ji,jj) = ztabramp(ji,jj) * ssmask(ji,jj) 288 fspf(ji,jj) = 0.25_wp * ( ztabramp(ji ,jj ) + ztabramp(ji ,jj+1) & 289 & +ztabramp(ji+1,jj+1) + ztabramp(ji+1,jj ) ) & 290 & * ssvmask(ji,jj) * ssvmask(ji,jj+1) 291 END DO 292 END DO 293 CALL lbc_lnk( 'agrif_Sponge', fspt, 'T', 1. ) ! Lateral boundary conditions 294 CALL lbc_lnk( 'agrif_Sponge', fspf, 'F', 1. ) 185 295 186 296 spongedoneU = .TRUE. 187 297 ENDIF 298 299 #if defined key_vertical 300 ! Remove vertical interpolation where not needed: 301 DO jj = 2, jpjm1 302 DO ji = 2, jpim1 303 IF ((fspu(ji-1,jj)==0._wp).AND.(fspu(ji,jj)==0._wp).AND. & 304 & (fspv(ji,jj-1)==0._wp).AND.(fspv(ji,jj)==0._wp)) mbkt_parent(ji,jj) = 0 305 ! 306 IF ((fspt(ji+1,jj)==0._wp).AND.(fspt(ji,jj)==0._wp).AND. & 307 & (fspf(ji,jj-1)==0._wp).AND.(fspf(ji,jj)==0._wp)) mbku_parent(ji,jj) = 0 308 ! 309 IF ((fspt(ji,jj+1)==0._wp).AND.(fspt(ji,jj)==0._wp).AND. & 310 & (fspf(ji-1,jj)==0._wp).AND.(fspf(ji,jj)==0._wp)) mbkv_parent(ji,jj) = 0 311 ! 312 IF ( ssmask(ji,jj) == 0._wp) mbkt_parent(ji,jj) = 0 313 IF (ssumask(ji,jj) == 0._wp) mbku_parent(ji,jj) = 0 314 IF (ssvmask(ji,jj) == 0._wp) mbkv_parent(ji,jj) = 0 315 END DO 316 END DO 317 ! 318 ztabramp(:,:) = REAL( mbkt_parent(:,:), wp ) ; CALL lbc_lnk( 'Agrif_Sponge', ztabramp, 'T', 1. ) 319 mbkt_parent(:,:) = NINT( ztabramp(:,:) ) 320 ztabramp(:,:) = REAL( mbku_parent(:,:), wp ) ; CALL lbc_lnk( 'Agrif_Sponge', ztabramp, 'U', 1. ) 321 mbku_parent(:,:) = NINT( ztabramp(:,:) ) 322 ztabramp(:,:) = REAL( mbkv_parent(:,:), wp ) ; CALL lbc_lnk( 'Agrif_Sponge', ztabramp, 'V', 1. ) 323 mbkv_parent(:,:) = NINT( ztabramp(:,:) ) 324 #endif 188 325 ! 189 326 #endif … … 201 338 INTEGER :: ji, jj, jk, jn ! dummy loop indices 202 339 INTEGER :: iku, ikv 203 REAL(wp) :: ztsa, zabe1, zabe2, zbtr 340 REAL(wp) :: ztsa, zabe1, zabe2, zbtr, zhtot, ztrelax 204 341 REAL(wp), DIMENSION(i1:i2,j1:j2,jpk) :: ztu, ztv 205 342 REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tsbdiff … … 210 347 REAL(wp), DIMENSION(1:jpk) :: h_out 211 348 INTEGER :: N_in, N_out 212 REAL(wp) :: h_diff213 349 !!---------------------------------------------------------------------- 214 350 ! … … 225 361 226 362 # 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 363 ! Interpolate thicknesses 364 ! Warning: these are masked, hence extrapolated prior interpolation. 365 DO jk=k1,k2 366 DO jj=j1,j2 367 DO ji=i1,i2 368 tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_b(ji,jj,jk) 369 END DO 370 END DO 371 END DO 372 373 ! Extrapolate thicknesses in partial bottom cells: 374 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 375 IF (ln_zps) THEN 376 DO jj=j1,j2 377 DO ji=i1,i2 378 jk = mbkt(ji,jj) 379 tabres(ji,jj,jk,jpts+1) = 0._wp 380 END DO 381 END DO 382 END IF 383 384 ! Save ssh at last level: 385 IF (.NOT.ln_linssh) THEN 386 tabres(i1:i2,j1:j2,k2,jpts+1) = sshb(i1:i2,j1:j2)*tmask(i1:i2,j1:j2,1) 387 ELSE 388 tabres(i1:i2,j1:j2,k2,jpts+1) = 0._wp 389 END IF 234 390 # endif 235 391 … … 237 393 ! 238 394 # if defined key_vertical 239 tabres_child(:,:,:,:) = 0. 395 396 IF (ln_linssh) tabres(i1:i2,j1:j2,k2,n2) = 0._wp 397 240 398 DO jj=j1,j2 241 399 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 400 tabres_child(ji,jj,:,:) = 0._wp 401 N_in = mbkt_parent(ji,jj) 402 zhtot = 0._wp 403 DO jk=1,N_in !k2 = jpk of parent grid 404 IF (jk==N_in) THEN 405 h_in(jk) = ht0_parent(ji,jj) + tabres(ji,jj,k2,n2) - zhtot 406 ELSE 407 h_in(jk) = tabres(ji,jj,jk,n2) 408 ENDIF 409 zhtot = zhtot + h_in(jk) 246 410 tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1) 247 h_in(N_in) = tabres(ji,jj,jk,n2)248 411 END DO 249 412 N_out = 0 … … 251 414 IF (tmask(ji,jj,jk) == 0) EXIT 252 415 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 above416 h_out(jk) = e3t_b(ji,jj,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above 254 417 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 418 419 ! Account for small differences in free-surface 420 IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN 421 h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) ) 422 ELSE 423 h_in(1) = h_in(1) - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) ) 424 ENDIF 425 IF (N_in*N_out > 0) THEN 426 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 427 ENDIF 262 428 ENDDO … … 268 434 DO jk=1,jpkm1 269 435 # if defined key_vertical 270 tsbdiff(ji,jj,jk,1:jpts) = tsb(ji,jj,jk,1:jpts) - tabres_child(ji,jj,jk,1:jpts)436 tsbdiff(ji,jj,jk,1:jpts) = (tsb(ji,jj,jk,1:jpts) - tabres_child(ji,jj,jk,1:jpts)) * tmask(ji,jj,jk) 271 437 # else 272 tsbdiff(ji,jj,jk,1:jpts) = tsb(ji,jj,jk,1:jpts) - tabres(ji,jj,jk,1:jpts)438 tsbdiff(ji,jj,jk,1:jpts) = (tsb(ji,jj,jk,1:jpts) - tabres(ji,jj,jk,1:jpts))*tmask(ji,jj,jk) 273 439 # endif 274 440 ENDDO 275 441 ENDDO 276 442 ENDDO 443 444 !* set relaxation time scale 445 IF( neuler == 0 .AND. lk_agrif_fstep ) THEN ; ztrelax = rn_trelax_tra / ( rdt ) 446 ELSE ; ztrelax = rn_trelax_tra / (2._wp * rdt ) 447 ENDIF 277 448 278 449 DO jn = 1, jpts … … 281 452 DO jj = j1,j2 282 453 DO ji = i1,i2-1 283 zabe1 = fsaht_spu(ji,jj) * umask(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk)454 zabe1 = rn_sponge_tra * fspu(ji,jj) * umask(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk) 284 455 ztu(ji,jj,jk) = zabe1 * ( tsbdiff(ji+1,jj ,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 285 456 END DO … … 288 459 DO ji = i1,i2 289 460 DO jj = j1,j2-1 290 zabe2 = fsaht_spv(ji,jj) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk)461 zabe2 = rn_sponge_tra * fspv(ji,jj) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk) 291 462 ztv(ji,jj,jk) = zabe2 * ( tsbdiff(ji ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) ) 292 463 END DO … … 312 483 zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) 313 484 ! horizontal diffusive trends 314 ztsa = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) 485 ztsa = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) & 486 & - ztrelax * fspt(ji,jj) * tsbdiff(ji,jj,jk,jn) 315 487 ! add it to the general tracer trends 316 488 tsa(ji,jj,jk,jn) = tsa(ji,jj,jk,jn) + ztsa … … 339 511 340 512 ! sponge parameters 341 REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, h_diff513 REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, zhtot, ztrelax 342 514 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ubdiff 343 515 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff … … 346 518 REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 347 519 REAL(wp), DIMENSION(1:jpk) :: h_out 348 INTEGER ::N_in, N_out520 INTEGER ::N_in, N_out 349 521 !!--------------------------------------------- 350 522 ! 351 523 IF( before ) THEN 352 DO jk= 1,jpkm1524 DO jk=k1,k2 353 525 DO jj=j1,j2 354 526 DO ji=i1,i2 355 527 tabres(ji,jj,jk,m1) = ub(ji,jj,jk) 356 528 # if defined key_vertical 357 tabres(ji,jj,jk,m2) = e3u_ n(ji,jj,jk)*umask(ji,jj,jk)529 tabres(ji,jj,jk,m2) = e3u_b(ji,jj,jk)*umask(ji,jj,jk) 358 530 # endif 359 531 END DO … … 361 533 END DO 362 534 535 # if defined key_vertical 536 ! Extrapolate thicknesses in partial bottom cells: 537 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 538 IF (ln_zps) THEN 539 DO jj=j1,j2 540 DO ji=i1,i2 541 jk = mbku(ji,jj) 542 tabres(ji,jj,jk,m2) = 0._wp 543 END DO 544 END DO 545 END IF 546 ! Save ssh at last level: 547 tabres(i1:i2,j1:j2,k2,m2) = 0._wp 548 IF (.NOT.ln_linssh) THEN 549 ! This vertical sum below should be replaced by the sea-level at U-points (optimization): 550 DO jk=1,jpk 551 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) 552 END DO 553 tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hu_0(i1:i2,j1:j2) 554 END IF 555 # endif 556 363 557 ELSE 364 558 365 559 # if defined key_vertical 366 tabres_child(:,:,:) = 0._wp 560 IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp 561 367 562 DO jj=j1,j2 368 563 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 564 tabres_child(ji,jj,:) = 0._wp 565 N_in = mbku_parent(ji,jj) 566 zhtot = 0._wp 567 DO jk=1,N_in 568 IF (jk==N_in) THEN 569 h_in(jk) = hu0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot 570 ELSE 571 h_in(jk) = tabres(ji,jj,jk,m2) 572 ENDIF 573 zhtot = zhtot + h_in(jk) 373 574 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 575 ENDDO 576 ! 577 N_out = 0 578 DO jk=1,jpk 579 IF (umask(ji,jj,jk) == 0) EXIT 580 N_out = N_out + 1 581 h_out(N_out) = e3u_b(ji,jj,jk) 582 ENDDO 583 584 ! Account for small differences in free-surface 585 IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN 586 h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) ) 587 ELSE 588 h_in(1) = h_in(1) - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) ) 589 ENDIF 590 591 IF (N_in * N_out > 0) THEN 592 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) 593 ENDIF 402 594 ENDDO 403 595 ENDDO … … 407 599 ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:) 408 600 #endif 601 !* set relaxation time scale 602 IF( neuler == 0 .AND. lk_agrif_fstep ) THEN ; ztrelax = rn_trelax_dyn / ( rdt ) 603 ELSE ; ztrelax = rn_trelax_dyn / (2._wp * rdt ) 604 ENDIF 409 605 ! 410 606 DO jk = 1, jpkm1 ! Horizontal slab … … 416 612 DO jj = j1,j2 417 613 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) ) * zbtr614 zbtr = r1_e1e2t(ji,jj) / e3t_b(ji,jj,jk) * rn_sponge_dyn * fspt(ji,jj) 615 hdivdiff(ji,jj,jk) = ( e2u(ji ,jj)*e3u_b(ji ,jj,jk) * ubdiff(ji ,jj,jk) & 616 & -e2u(ji-1,jj)*e3u_b(ji-1,jj,jk) * ubdiff(ji-1,jj,jk) ) * zbtr 421 617 END DO 422 618 END DO … … 424 620 DO jj = j1,j2-1 425 621 DO ji = i1,i2 ! vector opt. 426 zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * fsahm_spf(ji,jj)622 zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * rn_sponge_dyn * fspf(ji,jj) 427 623 rotdiff(ji,jj,jk) = ( -e1u(ji,jj+1) * ubdiff(ji,jj+1,jk) & 428 624 & +e1u(ji,jj ) * ubdiff(ji,jj ,jk) ) * fmask(ji,jj,jk) * zbtr … … 440 636 ! horizontal diffusive trends 441 637 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) 638 & + ( hdivdiff(ji+1,jj,jk) - ze1v ) * r1_e1u(ji,jj) & 639 & - ztrelax * fspu(ji,jj) * ubdiff(ji,jj,jk) 443 640 444 641 ! add it to the general momentum trends 445 ua(ji,jj,jk) = ua(ji,jj,jk) + zua 446 642 ua(ji,jj,jk) = ua(ji,jj,jk) + zua 447 643 END DO 448 644 ENDIF … … 492 688 ! 493 689 INTEGER :: ji, jj, jk, imax 494 REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, h_diff690 REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, zhtot, ztrelax 495 691 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: vbdiff 496 692 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff … … 503 699 504 700 IF( before ) THEN 505 DO jk= 1,jpkm1701 DO jk=k1,k2 506 702 DO jj=j1,j2 507 703 DO ji=i1,i2 508 704 tabres(ji,jj,jk,m1) = vb(ji,jj,jk) 509 705 # if defined key_vertical 510 tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v_ n(ji,jj,jk)706 tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v_b(ji,jj,jk) 511 707 # endif 512 708 END DO 513 709 END DO 514 710 END DO 711 712 # if defined key_vertical 713 ! Extrapolate thicknesses in partial bottom cells: 714 ! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on 715 IF (ln_zps) THEN 716 DO jj=j1,j2 717 DO ji=i1,i2 718 jk = mbkv(ji,jj) 719 tabres(ji,jj,jk,m2) = 0._wp 720 END DO 721 END DO 722 END IF 723 ! Save ssh at last level: 724 tabres(i1:i2,j1:j2,k2,m2) = 0._wp 725 IF (.NOT.ln_linssh) THEN 726 ! This vertical sum below should be replaced by the sea-level at V-points (optimization): 727 DO jk=1,jpk 728 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) 729 END DO 730 tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hv_0(i1:i2,j1:j2) 731 END IF 732 # endif 733 515 734 ELSE 516 735 517 736 # if defined key_vertical 518 tabres_child(:,:,:) = 0._wp737 IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp 519 738 DO jj=j1,j2 520 739 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 740 tabres_child(ji,jj,:) = 0._wp 741 N_in = mbkv_parent(ji,jj) 742 zhtot = 0._wp 743 DO jk=1,N_in 744 IF (jk==N_in) THEN 745 h_in(jk) = hv0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot 746 ELSE 747 h_in(jk) = tabres(ji,jj,jk,m2) 748 ENDIF 749 zhtot = zhtot + h_in(jk) 525 750 tabin(jk) = tabres(ji,jj,jk,m1) 526 h_in(N_in) = tabres(ji,jj,jk,m2) 527 ENDDO 751 ENDDO 752 ! 753 N_out = 0 754 DO jk=1,jpk 755 IF (vmask(ji,jj,jk) == 0) EXIT 756 N_out = N_out + 1 757 h_out(N_out) = e3v_b(ji,jj,jk) 758 ENDDO 759 760 ! Account for small differences in free-surface 761 IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN 762 h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) ) 763 ELSE 764 h_in(1) = h_in(1) - ( sum(h_in(1:N_in))-sum(h_out(1:N_out)) ) 765 ENDIF 528 766 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) 767 IF (N_in * N_out > 0) THEN 768 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) 769 ENDIF 548 770 ENDDO 549 771 ENDDO … … 553 775 vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:) 554 776 # endif 777 !* set relaxation time scale 778 IF( neuler == 0 .AND. lk_agrif_fstep ) THEN ; ztrelax = rn_trelax_dyn / ( rdt ) 779 ELSE ; ztrelax = rn_trelax_dyn / (2._wp * rdt ) 780 ENDIF 555 781 ! 556 782 DO jk = 1, jpkm1 ! Horizontal slab … … 562 788 DO jj = j1+1,j2 563 789 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) ) * zbtr790 zbtr = r1_e1e2t(ji,jj) / e3t_b(ji,jj,jk) * rn_sponge_dyn * fspt(ji,jj) 791 hdivdiff(ji,jj,jk) = ( e1v(ji,jj ) * e3v_b(ji,jj ,jk) * vbdiff(ji,jj ,jk) & 792 & -e1v(ji,jj-1) * e3v_b(ji,jj-1,jk) * vbdiff(ji,jj-1,jk) ) * zbtr 567 793 END DO 568 794 END DO 569 795 DO jj = j1,j2 570 796 DO ji = i1,i2-1 ! vector opt. 571 zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * fsahm_spf(ji,jj)797 zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * rn_sponge_dyn * fspf(ji,jj) 572 798 rotdiff(ji,jj,jk) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj,jk) & 573 799 & -e2v(ji ,jj) * vbdiff(ji ,jj,jk) ) * fmask(ji,jj,jk) * zbtr … … 602 828 va(ji,jj,jk) = va(ji,jj,jk) & 603 829 & + ( 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) 830 & + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji ,jj,jk) ) * r1_e2v(ji,jj) & 831 & - ztrelax * fspv(ji,jj) * vbdiff(ji,jj,jk) 605 832 END DO 606 833 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
