Changeset 7971 for branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90
- Timestamp:
- 2017-04-26T12:03:26+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r7963_nemo_v3_6_AGRIF-3_AGRIFVVL/NEMOGCM/NEMO/NST_SRC/agrif_opa_update.F90
r6204 r7971 13 13 USE dynspg_oce 14 14 USE zdf_oce ! vertical physics: ocean variables 15 USE domvvl ! Need interpolation routines 15 16 16 17 IMPLICIT NONE 17 18 PRIVATE 18 19 19 PUBLIC Agrif_Update_Tra, Agrif_Update_Dyn,Update_Scales 20 PUBLIC Agrif_Update_Tra, Agrif_Update_Dyn, Update_Scales, Agrif_Update_vvl 21 20 22 # if defined key_zdftke 21 23 PUBLIC Agrif_Update_Tke … … 177 179 # endif /* key_zdftke */ 178 180 181 RECURSIVE SUBROUTINE Agrif_Update_vvl( ) 182 !!--------------------------------------------- 183 !! *** ROUTINE Agrif_Update_vvl *** 184 !!--------------------------------------------- 185 ! 186 IF (Agrif_Root()) RETURN 187 ! 188 #if defined TWO_WAY 189 ! 190 IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Update e3 from grid Number',Agrif_Fixed(), 'Step', Agrif_Nb_Step() 191 ! 192 Agrif_UseSpecialValueInUpdate = .FALSE. 193 Agrif_SpecialValueFineGrid = 0. 194 ! 195 # if ! defined DECAL_FEEDBACK 196 CALL Agrif_Update_Variable(e3t_id, procname=updatee3t) 197 # else 198 CALL Agrif_Update_Variable(e3t_id, locupdate=(/1,0/), procname=updatee3t) 199 # endif 200 ! 201 Agrif_UseSpecialValueInUpdate = .FALSE. 202 ! 203 CALL Agrif_ChildGrid_To_ParentGrid() 204 CALL dom_vvl_update_UVF 205 CALL Agrif_ParentGrid_To_ChildGrid() 206 ! 207 IF ( lk_agrif_doupd ) THEN ! Initialisation: do recursive update: 208 CALL Agrif_ChildGrid_To_ParentGrid() 209 CALL Agrif_Update_vvl() 210 CALL Agrif_ParentGrid_To_ChildGrid() 211 ENDIF 212 ! 213 #endif 214 ! 215 END SUBROUTINE Agrif_Update_vvl 216 217 SUBROUTINE dom_vvl_update_UVF 218 !!--------------------------------------------- 219 !! *** ROUTINE dom_vvl_update_UVF *** 220 !!--------------------------------------------- 221 # include "domzgr_substitute.h90" 222 !! 223 INTEGER :: jk 224 REAL(wp):: zcoef 225 !!--------------------------------------------- 226 227 IF (lwp.AND.lk_agrif_debug) Write(*,*) 'Finalize e3 on grid Number', & 228 & Agrif_Fixed(), 'Step', Agrif_Nb_Step() 229 230 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' ) 231 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' ) 232 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' ) 233 234 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' ) 235 CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' ) 236 237 ! Update total depths: 238 hu(:,:) = 0._wp ! Ocean depth at U-points 239 hv(:,:) = 0._wp ! Ocean depth at V-points 240 DO jk = 1, jpkm1 241 hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk) 242 hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk) 243 END DO 244 ! ! Inverse of the local depth 245 hur(:,:) = 1._wp / ( hu(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:) 246 hvr(:,:) = 1._wp / ( hv(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:) 247 248 #if ! defined key_dynspg_ts 249 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 250 #else 251 IF ((.NOT.(lk_agrif_fstep.AND.(neuler==0))).AND.(.NOT.ln_bt_fw)) THEN 252 #endif 253 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' ) 254 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' ) 255 256 CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' ) 257 CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' ) 258 259 hu_b(:,:) = 0._wp ! Ocean depth at U-points 260 hv_b(:,:) = 0._wp ! Ocean depth at V-points 261 DO jk = 1, jpkm1 262 hu_b(:,:) = hu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk) 263 hv_b(:,:) = hv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 264 END DO 265 ! ! Inverse of the local depth 266 hur_b(:,:) = 1._wp / ( hu_b(:,:) + 1._wp - umask_i(:,:) ) * umask_i(:,:) 267 hvr_b(:,:) = 1._wp / ( hv_b(:,:) + 1._wp - vmask_i(:,:) ) * vmask_i(:,:) 268 ENDIF 269 270 ! 271 END SUBROUTINE dom_vvl_update_UVF 272 179 273 SUBROUTINE updateTS( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before ) 180 274 !!--------------------------------------------- 181 275 !! *** ROUTINE updateT *** 182 276 !!--------------------------------------------- 183 # include "domzgr_substitute.h90"184 277 INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2 185 278 REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres … … 188 281 INTEGER :: ji,jj,jk,jn 189 282 !!--------------------------------------------- 283 ! 190 284 ! 191 285 IF (before) THEN … … 194 288 DO jj=j1,j2 195 289 DO ji=i1,i2 196 tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) 290 !> jc tmp 291 tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) * fse3t_n(ji,jj,jk) / e3t_0(ji,jj,jk) 292 ! tabres(ji,jj,jk,jn) = tsn(ji,jj,jk,jn) * fse3t_n(ji,jj,jk) 293 !< jc tmp 197 294 END DO 198 295 END DO … … 200 297 END DO 201 298 ELSE 299 !> jc tmp 300 DO jn = n1,n2 301 tabres(i1:i2,j1:j2,k1:k2,jn) = tabres(i1:i2,j1:j2,k1:k2,jn) * e3t_0(i1:i2,j1:j2,k1:k2) 302 ENDDO 303 !< jc tmp 304 202 305 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 203 306 ! Add asselin part … … 209 312 tsb(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) & 210 313 & + atfp * ( tabres(ji,jj,jk,jn) & 211 & - tsn(ji,jj,jk,jn) ) * tmask(ji,jj,jk) 314 & - tsn(ji,jj,jk,jn)*fse3t_a(ji,jj,jk) ) & 315 & * tmask(ji,jj,jk) / fse3t_b(ji,jj,jk) 212 316 ENDIF 213 317 ENDDO … … 221 325 DO ji=i1,i2 222 326 IF( tabres(ji,jj,jk,jn) .NE. 0. ) THEN 223 tsn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) * tmask(ji,jj,jk) 327 #if defined key_vvl 328 tsn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) / fse3t_n(ji,jj,jk) 329 #else 330 tsn(ji,jj,jk,jn) = tabres(ji,jj,jk,jn) 331 #endif 224 332 END IF 225 333 END DO … … 260 368 DO jj=j1,j2 261 369 DO ji=i1,i2 262 tabres(ji,jj,jk) = tabres(ji,jj,jk) / e2u(ji,jj) / fse3u_n(ji,jj,jk)370 tabres(ji,jj,jk) = tabres(ji,jj,jk) / e2u(ji,jj) 263 371 ! 264 372 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 265 373 ub(ji,jj,jk) = ub(ji,jj,jk) & 266 & + atfp * ( tabres(ji,jj,jk) - un(ji,jj,jk) ) * umask(ji,jj,jk) 374 & + atfp * ( tabres(ji,jj,jk) & 375 & - un(ji,jj,jk)*fse3u_a(ji,jj,jk) ) & 376 & * umask(ji,jj,jk) / fse3u_b(ji,jj,jk) 267 377 ENDIF 268 378 ! 269 un(ji,jj,jk) = tabres(ji,jj,jk) * umask(ji,jj,jk) 379 un(ji,jj,jk) = tabres(ji,jj,jk) * umask(ji,jj,jk) / fse3u_n(ji,jj,jk) 270 380 END DO 271 381 END DO … … 304 414 DO jj=j1,j2 305 415 DO ji=i1,i2 306 tabres(ji,jj,jk) = tabres(ji,jj,jk) / e1v(ji,jj) / fse3v_n(ji,jj,jk)416 tabres(ji,jj,jk) = tabres(ji,jj,jk) / e1v(ji,jj) 307 417 ! 308 418 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 309 419 vb(ji,jj,jk) = vb(ji,jj,jk) & 310 & + atfp * ( tabres(ji,jj,jk) - vn(ji,jj,jk) ) * vmask(ji,jj,jk) 420 & + atfp * ( tabres(ji,jj,jk) & 421 & - vn(ji,jj,jk)*fse3v_a(ji,jj,jk) ) & 422 & * vmask(ji,jj,jk) / fse3v_b(ji,jj,jk) 311 423 ENDIF 312 424 ! 313 vn(ji,jj,jk) = tabres(ji,jj,jk) * vmask(ji,jj,jk) 425 vn(ji,jj,jk) = tabres(ji,jj,jk) * vmask(ji,jj,jk) / fse3v_n(ji,jj,jk) 314 426 END DO 315 427 END DO … … 345 457 DO jj=j1,j2 346 458 DO ji=i1,i2 347 tabres(ji,jj) = tabres(ji,jj) * hur(ji,jj)/ e2u(ji,jj)459 tabres(ji,jj) = tabres(ji,jj) / e2u(ji,jj) 348 460 ! 349 461 ! Update "now" 3d velocities: … … 352 464 spgu(ji,jj) = spgu(ji,jj) + fse3u_n(ji,jj,jk) * un(ji,jj,jk) 353 465 END DO 354 spgu(ji,jj) = spgu(ji,jj) * hur(ji,jj)355 466 ! 356 zcorr = tabres(ji,jj) - spgu(ji,jj)467 zcorr = (tabres(ji,jj) - spgu(ji,jj)) * hur(ji,jj) 357 468 DO jk=1,jpkm1 358 469 un(ji,jj,jk) = un(ji,jj,jk) + zcorr * umask(ji,jj,jk) … … 360 471 ! 361 472 ! Update barotropic velocities: 362 #if defined key_dynspg_ts 363 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 364 zcorr = tabres(ji,jj) - un_b(ji,jj) 473 #if ! defined key_dynspg_ts 474 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 475 #else 476 IF ((.NOT.(lk_agrif_fstep.AND.(neuler==0))).AND.(.NOT.ln_bt_fw)) THEN 477 #endif 478 zcorr = (tabres(ji,jj) - un_b(ji,jj) * hu_a(ji,jj)) * hur_b(ji,jj) 365 479 ub_b(ji,jj) = ub_b(ji,jj) + atfp * zcorr * umask(ji,jj,1) 366 END IF 367 #endif 368 un_b(ji,jj) = tabres(ji,jj) * umask(ji,jj,1) 480 END IF 481 un_b(ji,jj) = tabres(ji,jj) * hur(ji,jj) * umask(ji,jj,1) 369 482 ! 370 483 ! Correct "before" velocities to hold correct bt component: … … 373 486 spgu(ji,jj) = spgu(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk) 374 487 END DO 375 spgu(ji,jj) = spgu(ji,jj) * hur_b(ji,jj)376 488 ! 377 zcorr = ub_b(ji,jj) - spgu(ji,jj) 489 zcorr = ub_b(ji,jj) - spgu(ji,jj) * hur_b(ji,jj) 378 490 DO jk=1,jpkm1 379 491 ub(ji,jj,jk) = ub(ji,jj,jk) + zcorr * umask(ji,jj,jk) … … 410 522 DO jj=j1,j2 411 523 DO ji=i1,i2 412 tabres(ji,jj) = tabres(ji,jj) * hvr(ji,jj)/ e1v(ji,jj)524 tabres(ji,jj) = tabres(ji,jj) / e1v(ji,jj) 413 525 ! 414 526 ! Update "now" 3d velocities: … … 417 529 spgv(ji,jj) = spgv(ji,jj) + fse3v_n(ji,jj,jk) * vn(ji,jj,jk) 418 530 END DO 419 spgv(ji,jj) = spgv(ji,jj) * hvr(ji,jj)420 531 ! 421 zcorr = tabres(ji,jj) - spgv(ji,jj)532 zcorr = (tabres(ji,jj) - spgv(ji,jj)) * hvr(ji,jj) 422 533 DO jk=1,jpkm1 423 534 vn(ji,jj,jk) = vn(ji,jj,jk) + zcorr * vmask(ji,jj,jk) … … 425 536 ! 426 537 ! Update barotropic velocities: 427 #if defined key_dynspg_ts 428 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN ! Add asselin part 429 zcorr = tabres(ji,jj) - vn_b(ji,jj) 538 #if ! defined key_dynspg_ts 539 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 540 #else 541 IF ((.NOT.(lk_agrif_fstep.AND.(neuler==0))).AND.(.NOT.ln_bt_fw)) THEN 542 #endif 543 zcorr = (tabres(ji,jj) - vn_b(ji,jj)*hv_a(ji,jj)) * hvr_b(ji,jj) 430 544 vb_b(ji,jj) = vb_b(ji,jj) + atfp * zcorr * vmask(ji,jj,1) 431 END IF 432 #endif 433 vn_b(ji,jj) = tabres(ji,jj) * vmask(ji,jj,1) 545 END IF 546 vn_b(ji,jj) = tabres(ji,jj) * hvr(ji,jj) * vmask(ji,jj,1) 434 547 ! 435 548 ! Correct "before" velocities to hold correct bt component: … … 438 551 spgv(ji,jj) = spgv(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk) 439 552 END DO 440 spgv(ji,jj) = spgv(ji,jj) * hvr_b(ji,jj)441 553 ! 442 zcorr = vb_b(ji,jj) - spgv(ji,jj) 554 zcorr = vb_b(ji,jj) - spgv(ji,jj) * hvr_b(ji,jj) 443 555 DO jk=1,jpkm1 444 556 vb(ji,jj,jk) = vb(ji,jj,jk) + zcorr * vmask(ji,jj,jk) … … 471 583 ELSE 472 584 #if ! defined key_dynspg_ts 473 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 585 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 586 #else 587 IF ((.NOT.(lk_agrif_fstep.AND.(neuler==0))).AND.(.NOT.ln_bt_fw)) THEN 588 #endif 474 589 DO jj=j1,j2 475 590 DO ji=i1,i2 … … 479 594 END DO 480 595 ENDIF 481 #endif 596 ! 482 597 DO jj=j1,j2 483 598 DO ji=i1,i2 … … 655 770 # endif /* key_zdftke */ 656 771 772 SUBROUTINE updatee3t( ptab, i1, i2, j1, j2, k1, k2, before ) 773 !!--------------------------------------------- 774 !! *** ROUTINE updatee3t *** 775 !!--------------------------------------------- 776 INTEGER, INTENT(in) :: i1, i2, j1, j2, k1, k2 777 REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab 778 LOGICAL, INTENT(in) :: before 779 INTEGER :: ji,jj,jk 780 REAL(wp) :: zcoef 781 !!--------------------------------------------- 782 ! 783 IF (before) THEN 784 !> jc tmp: 785 ! ptab(i1:i2,j1:j2,k1:k2) = fse3t_n(i1:i2,j1:j2,k1:k2) 786 ptab(i1:i2,j1:j2,k1:k2) = fse3t_n(i1:i2,j1:j2,k1:k2) / e3t_0(i1:i2,j1:j2,k1:k2) 787 !< jc tmp: 788 ELSE 789 ! 790 ! 1) Updates at before time step: 791 ! ------------------------------- 792 ! 793 !> jc tmp: 794 ptab(i1:i2,j1:j2,k1:k2) = ptab(i1:i2,j1:j2,k1:k2) * e3t_0(i1:i2,j1:j2,k1:k2) 795 !< jc tmp: 796 797 #if ! defined key_dynspg_ts 798 IF (.NOT.(lk_agrif_fstep.AND.(neuler==0))) THEN 799 #else 800 IF ((.NOT.(lk_agrif_fstep.AND.(neuler==0))).AND.(.NOT.ln_bt_fw)) THEN 801 #endif 802 DO jk = 1, jpkm1 803 DO jj=j1,j2 804 DO ji=i1,i2 805 fse3t_b(ji,jj,jk) = fse3t_b(ji,jj,jk) & 806 & + atfp * ( ptab(ji,jj,jk) - fse3t_n(ji,jj,jk) ) 807 END DO 808 END DO 809 END DO 810 ! 811 fse3w_b (i1:i2,j1:j2,1) = e3w_0(i1:i2,j1:j2,1) + fse3t_b(i1:i2,j1:j2,1) - e3t_0(i1:i2,j1:j2,1) 812 fsdepw_b(i1:i2,j1:j2,1) = 0.0_wp 813 fsdept_b(i1:i2,j1:j2,1) = 0.5_wp * fse3w_b(i1:i2,j1:j2,1) 814 ! 815 DO jk = 2, jpkm1 816 DO jj = j1,j2 817 DO ji = i1,i2 818 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 819 fse3w_b(ji,jj,jk) = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * & 820 & ( fse3t_b(ji,jj,jk-1) - e3t_0(ji,jj,jk-1) ) & 821 & + 0.5_wp * tmask(ji,jj,jk) * & 822 & ( fse3t_b(ji,jj,jk ) - e3t_0(ji,jj,jk ) ) 823 fsdepw_b(ji,jj,jk) = fsdepw_b(ji,jj,jk-1) + fse3t_b(ji,jj,jk-1) 824 fsdept_b(ji,jj,jk) = zcoef * ( fsdepw_b(ji,jj,jk ) + 0.5 * fse3w_b(ji,jj,jk)) & 825 & + (1-zcoef) * ( fsdept_b(ji,jj,jk-1) + fse3w_b(ji,jj,jk)) 826 END DO 827 END DO 828 END DO 829 ENDIF 830 ! 831 ! 2) Updates at now time step: 832 ! ---------------------------- 833 ! 834 ! Update vertical scale factor at T-points: 835 fse3t_n(i1:i2,j1:j2,k1:k2) = ptab(i1:i2,j1:j2,k1:k2) 836 ! 837 ! Update total depth: 838 ht(i1:i2,j1:j2) = 0._wp 839 DO jk = 1, jpkm1 840 ht(i1:i2,j1:j2) = ht(i1:i2,j1:j2) + fse3t_n(i1:i2,j1:j2,jk) * tmask(i1:i2,j1:j2,jk) 841 END DO 842 ! 843 ! Update vertical scale factor at W-points and depths: 844 fse3w_n (i1:i2,j1:j2,1) = e3w_0(i1:i2,j1:j2,1) + fse3t_n(i1:i2,j1:j2,1) - e3t_0(i1:i2,j1:j2,1) 845 fsdept_n(i1:i2,j1:j2,1) = 0.5_wp * fse3w_n(i1:i2,j1:j2,1) 846 fsdepw_n(i1:i2,j1:j2,1) = 0.0_wp 847 fsde3w_n(i1:i2,j1:j2,1) = fsdept_n(i1:i2,j1:j2,1) - (ht(i1:i2,j1:j2)-ht_0(i1:i2,j1:j2)) ! Last term in the rhs is ssh 848 ! 849 DO jk = 2, jpkm1 850 DO jj = j1,j2 851 DO ji = i1,i2 852 zcoef = (tmask(ji,jj,jk) - wmask(ji,jj,jk)) 853 fse3w_n(ji,jj,jk) = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * tmask(ji,jj,jk) ) * ( fse3t_n(ji,jj,jk-1) - e3t_0(ji,jj,jk-1) ) & 854 & + 0.5_wp * tmask(ji,jj,jk) * ( fse3t_n(ji,jj,jk ) - e3t_0(ji,jj,jk ) ) 855 fsdepw_n(ji,jj,jk) = fsdepw_n(ji,jj,jk-1) + fse3t_n(ji,jj,jk-1) 856 fsdept_n(ji,jj,jk) = zcoef * ( fsdepw_n(ji,jj,jk ) + 0.5 * fse3w_n(ji,jj,jk)) & 857 & + (1-zcoef) * ( fsdept_n(ji,jj,jk-1) + fse3w_n(ji,jj,jk)) 858 fsde3w_n(ji,jj,jk) = fsdept_n(ji,jj,jk) - (ht(ji,jj)-ht_0(ji,jj)) ! Last term in the rhs is ssh 859 END DO 860 END DO 861 END DO 862 ! 863 ENDIF 864 ! 865 END SUBROUTINE updatee3t 866 657 867 #else 658 868 CONTAINS
Note: See TracChangeset
for help on using the changeset viewer.
