Changeset 8835
- Timestamp:
- 2017-11-28T15:16:43+01:00 (6 years ago)
- Location:
- branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC
- Files:
-
- 2 edited
-
agrif_opa_sponge.F90 (modified) (5 diffs)
-
agrif_user.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90
r7646 r8835 1 # define SPONGE && defineSPONGE_TOP1 #undef SPONGE && undef SPONGE_TOP 2 2 3 3 MODULE agrif_opa_sponge … … 59 59 REAL(wp) :: timecoeff 60 60 !!--------------------------------------------- 61 62 61 #if defined SPONGE 63 62 timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot() … … 203 202 INTEGER :: iku, ikv 204 203 REAL(wp) :: ztsa, zabe1, zabe2, zbtr 205 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ztu, ztv 206 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2) ::tsbdiff 204 REAL(wp), DIMENSION(i1:i2,j1:j2,jpk) :: ztu, ztv 205 REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tsbdiff 206 REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tabres_child 207 REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin 208 REAL(wp) :: h_in(k1:k2) 209 REAL(wp) :: h_out(1:jpk) 210 INTEGER :: N_in, N_out 211 REAL(wp) :: h_diff, zrhoxy 207 212 ! 208 213 IF( before ) THEN 209 tabres(i1:i2,j1:j2,k1:k2,n1:n2) = tsn(i1:i2,j1:j2,k1:k2,n1:n2) 214 DO jn = n1,n2-1 215 DO jk=k1,k2 216 DO jj=j1,j2 217 DO ji=i1,i2 218 tabres(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) 219 END DO 220 END DO 221 END DO 222 END DO 223 DO jk=k1,k2 224 DO jj=j1,j2 225 DO ji=i1,i2 226 tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) 227 END DO 228 END DO 229 END DO 210 230 ELSE 211 231 ! 212 tsbdiff(:,:,:,:) = tsb(i1:i2,j1:j2,:,:) - tabres(:,:,:,:) 232 #ifdef key_vertical 233 tabres_child(:,:,:,:) = 0. 234 zrhoxy = Agrif_rhox()*Agrif_rhoy() 235 DO jj=j1,j2 236 DO ji=i1,i2 237 N_in = 0 238 DO jk=k1,k2 !k2 = jpk of parent grid 239 IF (tabres(ji,jj,jk,n2) == 0) EXIT 240 N_in = N_in + 1 241 tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) 242 h_in(N_in) = tabres(ji,jj,jk,n2)/(e1e2t(ji,jj)*zrhoxy) 243 END DO 244 N_out = 0 245 DO jk=1,jpk ! jpk of child grid 246 IF (tmask(ji,jj,jk) == 0) EXIT ! TODO: Will not work with ISF. !This doesn't seem to work at the moment in GYRE but is OK in overflow model 247 N_out = N_out + 1 248 h_out(jk) = e3t_n(ji,jj,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above 249 ENDDO 250 IF (N_in > 0) THEN 251 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 252 tabres(ji,jj,k2,:) = tabres(ji,jj,k2-1,:) !what is this line for????? 253 DO jn=1,jpts 254 call reconstructandremap(tabin(1:N_in,jn),h_in,tabres_child(ji,jj,1:N_out,jn),h_out,N_in,N_out) 255 ENDDO 256 ENDIF 257 DO jk=1,jpk 258 tsbdiff(ji,jj,jk,n1:n2-1) = tsb(ji, jj,jk,n1:n2-1) - tabres_child(ji,jj,jk,n1:n2-1) 259 ENDDO 260 ENDDO 261 ENDDO 262 #else 263 do jk=k1,k2 264 do jj=j1,j2 265 do ji=i1,i2 266 ! This will be slow - Is there a way to speed it up and avoid divide by zero? 267 IF (tabres(ji,jj,jk,n2) .NE. 0) THEN 268 tsbdiff(ji,jj,jk,n1:n2-1) = tsb(ji,jj,jk,n1:n2) - tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) 269 ELSE 270 tsbdiff(ji,jj,jk,n1:n2-1) = 0._wp 271 ENDIF 272 enddo 273 enddo 274 enddo 275 #endif 213 276 DO jn = 1, jpts 214 277 DO jk = 1, jpkm1 … … 258 321 259 322 260 SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2, before)323 SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before) 261 324 !!--------------------------------------------- 262 325 !! *** ROUTINE interpun_sponge *** 263 326 !!--------------------------------------------- 264 INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2 265 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2 ), INTENT(inout) :: tabres327 INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2 328 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres 266 329 LOGICAL, INTENT(in) :: before 267 330 268 INTEGER :: ji,jj,jk 331 INTEGER :: ji,jj,jk,N_in,N_out 269 332 270 333 ! sponge parameters 271 REAL(wp) :: ze2u, ze1v, zua, zva, zbtr 272 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ubdiff 273 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: rotdiff, hdivdiff 334 REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, h_diff, zrhoy 335 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ubdiff 336 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff 337 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child 338 REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 339 REAL(wp), DIMENSION(1:jpk) :: h_out 274 340 INTEGER :: jmax 275 341 !!--------------------------------------------- 276 342 ! 277 343 IF( before ) THEN 278 tabres = un(i1:i2,j1:j2,:) 344 #ifdef key_vertical 345 DO jk=1,jpk 346 DO jj=j1,j2 347 DO ji=i1,i2 348 tabres(ji,jj,jk,m1) = (e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk)*umask(ji,jj,jk)) 349 tabres(ji,jj,jk,m2) = (umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk)) 350 END DO 351 END DO 352 END DO 353 354 #else 355 tabres(i1:i2,j1:j2,:,m1) = un(i1:i2,j1:j2,:) 356 #endif 279 357 ELSE 358 #ifdef key_vertical 359 tabres_child(:,:,:) = 0._wp 360 zrhoy = Agrif_rhoy() 361 DO jj=j1,j2 362 DO ji=i1,i2 363 N_in = 0 364 DO jk=k1,k2 365 IF (tabres(ji,jj,jk,m2) == 0) EXIT 366 N_in = N_in + 1 367 tabin(jk) = tabres(ji,jj,jk,m1)/tabres(ji,jj,jk,m2) 368 h_in(N_in) = tabres(ji,jj,jk,m2)/(e2u(ji,jj)*zrhoy) 369 ENDDO 370 ! 371 IF (N_in == 0) THEN 372 tabres_child(ji,jj,:) = 0. 373 CYCLE 374 ENDIF 375 376 N_out = 0 377 DO jk=1,jpk 378 if (umask(ji,jj,jk) == 0) EXIT 379 N_out = N_out + 1 380 h_out(N_out) = e3u_n(ji,jj,jk) 381 ENDDO 382 383 IF (N_out == 0) THEN 384 tabres_child(ji,jj,:) = 0. 385 CYCLE 386 ENDIF 387 388 IF (N_in * N_out > 0) THEN 389 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 390 if (h_diff < -1.e4) then 391 print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in)) 392 endif 393 ENDIF 394 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) 395 396 ENDDO 397 ENDDO 398 DO jk = 1, jpkm1 399 DO jj=j1,j2 400 ubdiff(i1:i2,jj,jk) = (ub(i1:i2,jj,jk) - tabres_child(i1:i2,jj,jk))*umask(i1:i2,jj,jk) 401 END DO 402 END DO 403 #else 280 404 ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres(:,:,:))*umask(i1:i2,j1:j2,:) 405 #endif 281 406 ! 282 407 DO jk = 1, jpkm1 ! Horizontal slab … … 356 481 357 482 358 SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2, before,nb,ndir)483 SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before,nb,ndir) 359 484 !!--------------------------------------------- 360 485 !! *** ROUTINE interpvn_sponge *** 361 486 !!--------------------------------------------- 362 INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2 363 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2 ), INTENT(inout) :: tabres487 INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2 488 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres 364 489 LOGICAL, INTENT(in) :: before 365 490 INTEGER, INTENT(in) :: nb , ndir 366 491 ! 367 INTEGER :: ji, jj, jk 368 REAL(wp) :: ze2u, ze1v, zua, zva, zbtr 369 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: vbdiff 370 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: rotdiff, hdivdiff 492 INTEGER :: ji, jj, jk, N_in, N_out 493 REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, h_diff, zrhox 494 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: vbdiff 495 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff 496 REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child 497 REAL(wp), DIMENSION(k1:k2) :: tabin, h_in 498 REAL(wp), DIMENSION(1:jpk) :: h_out 371 499 INTEGER :: imax 372 500 !!--------------------------------------------- 373 501 374 502 IF( before ) THEN 503 #ifdef key_vertical 504 DO jk=1,jpk 505 DO jj=j1,j2 506 DO ji=i1,i2 507 tabres(ji,jj,jk,m1) = (e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)*vmask(ji,jj,jk)) 508 tabres(ji,jj,jk,m2) = (vmask(ji,jj,jk) * e1v(ji,jj) * e3v_n(ji,jj,jk)) 509 END DO 510 END DO 511 END DO 512 #else 375 513 tabres = vn(i1:i2,j1:j2,:) 514 #endif 376 515 ELSE 516 #ifdef key_vertical 517 zrhox = Agrif_rhox() 518 tabres_child(:,:,:) = 0._wp 519 DO jj=j1,j2 520 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 525 tabin(jk) = tabres(ji,jj,jk,m1)/tabres(ji,jj,jk,m2) 526 h_in(N_in) = tabres(ji,jj,jk,m2)/(e1v(ji,jj)*zrhox) 527 ENDDO 528 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 (umask(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) 548 ENDDO 549 ENDDO 550 DO jk = 1, jpkm1 551 DO jj=j1,j2 552 vbdiff(i1:i2,jj,jk) = (vb(i1:i2,jj,jk) - tabres_child(i1:i2,jj,jk))*vmask(i1:i2,jj,jk) 553 END DO 554 END DO 555 #else 377 556 ! 378 557 vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres(:,:,:))*vmask(i1:i2,j1:j2,:) 558 #endif 379 559 ! 380 560 DO jk = 1, jpkm1 ! Horizontal slab -
branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_user.F90
r8135 r8835 347 347 !--------------------------------------------------------------------- 348 348 CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts+1/),tsn_id) 349 CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts /),tsn_sponge_id)349 CALL agrif_declare_variable((/2,2,0,0/),(/3,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,jpts+1/),tsn_sponge_id) 350 350 351 351 CALL agrif_declare_variable((/1,2,0,0/),(/2,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_interp_id) … … 353 353 CALL agrif_declare_variable((/1,2,0,0/),(/2,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_update_id) 354 354 CALL agrif_declare_variable((/2,1,0,0/),(/3,2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),vn_update_id) 355 CALL agrif_declare_variable((/1,2,0 /),(/2,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),un_sponge_id)356 CALL agrif_declare_variable((/2,1,0 /),(/3,2,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),vn_sponge_id)355 CALL agrif_declare_variable((/1,2,0,0/),(/2,3,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),un_sponge_id) 356 CALL agrif_declare_variable((/2,1,0,0/),(/3,2,0,0/),(/'x','y','N','N'/),(/1,1,1,1/),(/nlci,nlcj,jpk,2/),vn_sponge_id) 357 357 358 358 CALL agrif_declare_variable((/2,2,0/),(/3,3,0/),(/'x','y','N'/),(/1,1,1/),(/nlci,nlcj,jpk/),e3t_id)
Note: See TracChangeset
for help on using the changeset viewer.
