Changeset 8822 for branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90
- Timestamp:
- 2017-11-27T14:55:00+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90
r8690 r8822 206 206 DO jj=j1,j2 207 207 DO ji=i1,i2 208 tabres(ji,jj,jk,jn) = ( zrho_xy *tsn(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) &209 * tmask(ji,jj,jk) + (tmask(ji,jj,jk)-1)*999._wp208 tabres(ji,jj,jk,jn) = (tsn(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) ) & 209 * tmask(ji,jj,jk) * zrho_xy + (tmask(ji,jj,jk)-1)*999._wp 210 210 END DO 211 211 END DO … … 215 215 DO jj=j1,j2 216 216 DO ji=i1,i2 217 tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) &217 tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e1e2t(ji,jj) * e3t_n(ji,jj,jk) * zrho_xy & 218 218 + (tmask(ji,jj,jk)-1)*999._wp 219 219 END DO … … 240 240 N_in = 0 241 241 DO jk=k1,k2 !k2 = jpk of child grid 242 IF (tabres(ji,jj,jk,n2) > -900) EXIT242 IF (tabres(ji,jj,jk,n2) == 0 ) EXIT 243 243 N_in = N_in + 1 244 244 tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2) … … 247 247 N_out = 0 248 248 DO jk=1,jpk ! jpk of parent grid 249 IF (tmask(ji,jj,jk) ==0) EXIT ! TODO: Will not work with ISF249 IF (tmask(ji,jj,jk) < -900) EXIT ! TODO: Will not work with ISF 250 250 N_out = N_out + 1 251 251 h_out(N_out) = e3t_n(ji,jj,jk) !Parent grid scale factors. Could multiply by e1e2t here instead of division above 252 252 ENDDO 253 ! IF(ji.EQ.i1 .AND. jj.EQ.j1) print *,'1st parent point',sum(h_in(1:N_in)), sum(h_out(1:N_out))254 253 IF (N_in > 0) THEN 255 254 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 256 ! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly257 ! IF abs(h_diff > 1.e-8) THEN258 ! N_in = N_in + 1259 ! h_in(N_in) = h_diff260 ! tabin(N_in,:) = tabin(N_in-1,:)261 255 IF (h_diff < -1.e-4) THEN 262 print *,'CHECK YOUR bathy T points ...',ji,jj,h_diff,e1e2t(ji,jj),sum(h_in(1:N_in)),sum(h_out(1:N_out)), N_in, N_out 263 print *,h_in(1:N_in) 264 print *,h_out(1:N_out) 265 STOP 266 ! N_out = N_out + 1 267 ! h_out(N_out) = - h_diff 256 print *,'CHECK YOUR bathy T points ...',ji,jj,h_diff,e1e2t(ji,jj),sum(h_in(1:N_in)),sum(h_out(1:N_out)), N_in, N_out 257 print *,h_in(1:N_in) 258 print *,h_out(1:N_out) 259 STOP 268 260 ENDIF 269 261 DO jn=n1,n2-1 … … 276 268 tabres_child(:,:,:,:) = tabres(:,:,:,1:jpts) 277 269 #endif 278 279 270 ! WARNING : 280 271 ! tabres replaced by tabres_child in the following 281 272 ! k1 k2 replaced by 1 jpk 282 273 ! VERTICAL REFINEMENT END 283 284 274 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 285 275 ! Add asselin part … … 329 319 REAL(wp) :: h_out(1:jpk) 330 320 INTEGER :: N_in, N_out 331 REAL(wp) :: h_diff 321 REAL(wp) :: h_diff, excess, thick 332 322 REAL(wp) :: tabin(k1:k2) 333 323 ! VERTICAL REFINEMENT END … … 343 333 DO ji=i1,i2 344 334 tabres(ji,jj,jk,1) = zrhoy * e2u(ji,jj) * e3u_n(ji,jj,jk) * umask(ji,jj,jk) * un(ji,jj,jk) & 345 -(umask(ji,jj,jk)-1)*999._wp346 tabres(ji,jj,jk,2) = umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk) &347 -(umask(ji,jj,jk)-1)*999._wp335 + (umask(ji,jj,jk)-1)*999._wp 336 tabres(ji,jj,jk,2) = zrhoy * umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk) & 337 + (umask(ji,jj,jk)-1)*999._wp 348 338 END DO 349 339 END DO … … 357 347 tabres_child(:,:,:) = 0. 358 348 # if defined key_vertical 359 AGRIF_SpecialValue = -999._wp349 AGRIF_SpecialValue = 0._wp 360 350 DO jj=j1,j2 361 351 DO ji=i1,i2 362 352 N_in = 0 353 h_in(:) = 0._wp 354 tabin(:) = 0._wp 363 355 DO jk=k1,k2 !k2=jpk of child grid 364 IF (tabres(ji,jj,jk,2) >-900) EXIT356 IF( tabres(ji,jj,jk,2) < -900) EXIT 365 357 N_in = N_in + 1 366 358 tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) … … 375 367 IF (N_in * N_out > 0) THEN 376 368 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 377 ! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly 378 if (h_diff < -1.e-4) then 379 print *,'CHECK YOUR bathy U points ...',ji,jj,h_diff,e2u(ji,jj),sum(h_in(1:N_in)),sum(h_out(1:N_out)), N_in, N_out 380 stop 381 ! else ! Extends with 0 382 ! N_in = N_in + 1 383 ! tabin(N_in) = 0. 384 ! h_in(N_in) = h_diff 385 endif 369 IF (h_diff < -1.e-4) then 370 !Even if bathy at T points match it's possible for the U points to be deeper in the child grid. 371 !In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell. 372 excess = 0._wp 373 DO jk=N_in,1,-1 374 thick = MIN(-1*h_diff, h_in(jk)) 375 excess = excess + tabin(jk)*thick*e2u(ji,jj) 376 tabin(jk) = tabin(jk)*(1. - thick/h_in(jk)) 377 h_diff = h_diff + thick 378 IF ( h_diff == 0) THEN 379 N_in = jk 380 h_in(jk) = h_in(jk) - thick 381 EXIT 382 ENDIF 383 ENDDO 384 ENDIF 386 385 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) 387 ENDIF 386 tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e2u(ji,jj)*h_out(N_out)) 387 ENDIF 388 388 ENDDO 389 389 ENDDO 390 390 #else 391 391 DO jk=1,jpk … … 397 397 END DO 398 398 #endif 399 400 ! WARNING :401 ! tabres replaced by tabres_child in the following402 ! k1 k2 replaced by 1 jpk403 ! remove division by fs e3u (already done)404 ! VERTICAL REFINEMENT END405 399 406 400 DO jk=1,jpk 407 401 DO jj=j1,j2 408 402 DO ji=i1,i2 409 !Following line now replaced by division higher up in vertical refinement case410 ! tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e2u(ji,jj) / e3u_n(ji,jj,jk)411 !412 403 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 413 404 ub(ji,jj,jk) = ub(ji,jj,jk) & … … 439 430 REAL(wp) :: h_out(1:jpk) 440 431 INTEGER :: N_in, N_out 441 REAL(wp) :: h_diff 432 REAL(wp) :: h_diff, excess, thick 442 433 REAL(wp) :: tabin(k1:k2) 443 434 ! VERTICAL REFINEMENT END … … 447 438 zrhox = Agrif_Rhox() 448 439 #if defined key_vertical 440 AGRIF_SpecialValue = -999._wp 449 441 DO jk=k1,k2 450 442 DO jj=j1,j2 451 443 DO ji=i1,i2 452 tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk) 453 tabres(ji,jj,jk,2) = vmask(ji,jj,jk) * zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) 444 tabres(ji,jj,jk,1) = zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) * vmask(ji,jj,jk) * vn(ji,jj,jk) & 445 + (vmask(ji,jj,jk)-1)*999._wp 446 tabres(ji,jj,jk,2) = vmask(ji,jj,jk) * zrhox * e1v(ji,jj) * e3v_n(ji,jj,jk) & 447 + (vmask(ji,jj,jk)-1)*999._wp 454 448 END DO 455 449 END DO … … 466 460 ELSE 467 461 tabres_child(:,:,:) = 0. 468 ! VERTICAL REFINEMENT BEGIN469 462 #if defined key_vertical 463 AGRIF_SpecialValue = 0._wp 470 464 DO jj=j1,j2 471 465 DO ji=i1,i2 472 466 N_in = 0 473 467 DO jk=k1,k2 474 IF (tabres(ji,jj,jk,2) ==0) EXIT468 IF (tabres(ji,jj,jk,2) < -900) EXIT 475 469 N_in = N_in + 1 476 470 tabin(jk) = tabres(ji,jj,jk,1)/tabres(ji,jj,jk,2) … … 485 479 IF (N_in * N_out > 0) THEN 486 480 h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in)) 487 ! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly 488 if (h_diff < -1.e-4) then 489 print *,'CHECK YOUR BATHY ...' 490 stop 491 ! else ! Extends with 0 492 ! N_in = N_in + 1 493 ! tabin(N_in) = 0. 494 ! h_in(N_in) = h_diff 495 endif 481 IF (h_diff < -1.e-4) then 482 !Even if bathy at T points match it's possible for the U points to be deeper in the child grid. 483 !In this case we need to move transport from the child grid cells below bed of parent grid into the bottom cell. 484 excess = 0._wp 485 DO jk=N_in,1,-1 486 thick = MIN(-1*h_diff, h_in(jk)) 487 excess = excess + tabin(jk)*thick*e2u(ji,jj) 488 tabin(jk) = tabin(jk)*(1. - thick/h_in(jk)) 489 h_diff = h_diff + thick 490 IF ( h_diff == 0) THEN 491 N_in = jk 492 h_in(jk) = h_in(jk) - thick 493 EXIT 494 ENDIF 495 ENDDO 496 ENDIF 496 497 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) 498 tabres_child(ji,jj,N_out) = tabres_child(ji,jj,N_out) + excess/(e1v(ji,jj)*h_out(N_out)) 497 499 ENDIF 498 500 ENDDO … … 507 509 END DO 508 510 #endif 509 510 ! WARNING :511 ! tabres replaced by tabres_child in the following512 ! k1 k2 replaced by 1 jpk513 ! remove division by fs e3v (already done)514 ! VERTICAL REFINEMENT END515 511 516 512 DO jk=k1,jpk 517 513 DO jj=j1,j2 518 514 DO ji=i1,i2 519 !Following line now replaced by division higher up in vertical refinement case520 ! tabres(ji,jj,jk) = tabres(ji,jj,jk) * r1_e1v(ji,jj) / e3v_n(ji,jj,jk)521 515 ! 522 516 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part … … 761 755 END DO 762 756 ENDIF 763 !764 757 END SUBROUTINE updatevb2b 765 758
Note: See TracChangeset
for help on using the changeset viewer.
