Changeset 4194 for branches/2013
- Timestamp:
- 2013-11-14T12:18:20+01:00 (10 years ago)
- Location:
- branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN
- Files:
-
- 2 edited
-
dynspg_ts.F90 (modified) (13 diffs)
-
dynzdf_imp.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r4178 r4194 172 172 REAL(wp) :: zx1, zy1, zx2, zy2 ! - - 173 173 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 174 REAL(wp) :: za 1, za2, za3! - -174 REAL(wp) :: za0, za1, za2, za3 ! - - 175 175 REAL(wp) :: zu_spg, zu_cor, zu_sld, zu_asp ! - - 176 176 REAL(wp) :: zv_spg, zv_cor, zv_sld, zv_asp ! - - … … 182 182 REAL(wp), POINTER, DIMENSION(:,:) :: zhu_b, zhv_b 183 183 REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zsshp2_e 184 REAL(wp), POINTER, DIMENSION(:,:) :: zhust_e, zhvst_e 184 185 !!---------------------------------------------------------------------- 185 186 ! … … 191 192 CALL wrk_alloc( jpi, jpj, zhu_b, zhv_b ) 192 193 CALL wrk_alloc( jpi, jpj, zhup2_e, zhvp2_e, zsshp2_e ) 194 CALL wrk_alloc( jpi, jpj, zhust_e, zhvst_e ) 193 195 ! 194 196 … … 385 387 zcoef = -1._wp * z1_2dt_b 386 388 387 !--------> MERGED TO HERE388 389 IF(ln_bfrimp) THEN 389 390 ! ! Remove the bottom stress trend from 3-D sea surface level gradient … … 455 456 ! Initialize depths: 456 457 IF ( lk_vvl.AND.(.NOT.ln_bt_fw) ) THEN 457 ! hu_e (:,:) = hu_0(:,:) + sshu_b(:,:) 458 ! hv_e (:,:) = hv_0(:,:) + sshv_b(:,:) 459 ! hur_e (:,:) = umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) ) 460 ! hvr_e (:,:) = vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) ) 461 hu_e (:,:) = hu (:,:) 462 hv_e (:,:) = hv (:,:) 463 hur_e (:,:) = hur (:,:) 464 hvr_e (:,:) = hvr (:,:) 458 hur_e (:,:) = zhu_b (:,:) 459 hvr_e (:,:) = zhv_b (:,:) 460 hu_e(:,:) = umask(:,:,1) / ( zhu_b(:,:) + 1. - umask(:,:,1) ) 461 hv_e(:,:) = vmask(:,:,1) / ( zhv_b(:,:) + 1. - vmask(:,:,1) ) 465 462 ELSE 466 463 hu_e (:,:) = hu (:,:) … … 526 523 527 524 ! Extrapolate barotropic velocities at step jit+0.5: 528 ua_e(:,:) = za1 * zun_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)529 va_e(:,:) = za1 * zvn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)525 ua_e(:,:) = za1 * zun_e(:,:) + za2 * zub_e(:,:) + za3 * ubb_e(:,:) 526 va_e(:,:) = za1 * zvn_e(:,:) + za2 * zvb_e(:,:) + za3 * vbb_e(:,:) 530 527 531 528 IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only) 532 529 ! ! ------------------ 533 530 ! Extrapolate Sea Level at step jit+0.5: 534 zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:)531 zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * zsshb_e(:,:) + za3 * sshbb_e(:,:) 535 532 ! 536 533 DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points … … 554 551 ! considering fluxes below: 555 552 ! 556 zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5553 zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5 ! ACC INSTANT FAILURE IF THESE ARE USED (VERY LARGE NEGATIVE SALINITY AFTER 1 TS) 557 554 zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 555 ! zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:) ! fluxes at jn+0.5 ! ACC WORKS BETTER BUT FAILS AFTER O(100) TIMESTEPS 556 ! zwy(:,:) = e1v(:,:) * zvn_e(:,:) * hv_e(:,:) 558 557 DO jj = 2, jpjm1 559 558 DO ji = fs_2, fs_jpim1 ! vector opt. 560 zhdiv(ji,jj) = ( e2u(ji ,jj) * zun_e(ji ,jj) * hu_e(ji ,jj) & 561 & - e2u(ji-1,jj) * zun_e(ji-1,jj) * hu_e(ji-1,jj) & 562 & + e1v(ji,jj ) * zvn_e(ji,jj ) * hv_e(ji,jj ) & 563 & - e1v(ji,jj-1) * zvn_e(ji,jj-1) * hv_e(ji,jj-1) ) / ( e1t(ji,jj) * e2t(ji,jj) ) 559 zhdiv(ji,jj) = ( zwx(ji,jj) - zwx(ji-1,jj) & 560 & + zwy(ji,jj) - zwy(ji,jj-1) & 561 & ) / ( e1t(ji,jj) * e2t(ji,jj) ) 564 562 END DO 565 563 END DO … … 582 580 END DO 583 581 END DO 584 585 ! !* after barotropic velocities (vorticity scheme dependent) 586 ! ! --------------------------- 587 zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:) ! now_e transport 588 zwy(:,:) = e1v(:,:) * zvn_e(:,:) * hv_e(:,:) 582 #if defined key_bdy 583 ! Duplicate sea level across open boundaries (this is only cosmetic if lk_vvl=.false.) 584 IF (lk_bdy) CALL bdy_ssh( ssha_e ) 585 #endif 586 ! 587 ! Half-step back interpolation of SSH for surface pressure computation: 588 !---------------------------------------------------------------------- 589 IF ((jn==1).AND.ll_init) THEN 590 za0=1._wp ! Forward-backward 591 za1=0._wp 592 za2=0._wp 593 za3=0._wp 594 ELSEIF ((jn==2).AND.ll_init) THEN ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12 595 za0= 1.0833333333333_wp ! za0 = 1-gam-eps 596 za1=-0.1666666666666_wp ! za1 = gam 597 za2= 0.0833333333333_wp ! za2 = eps 598 za3= 0._wp 599 ELSE ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880 600 za0=0.614_wp ! za0 = 1/2 + gam + 2*eps 601 za1=0.285_wp ! za1 = 1/2 - 2*gam - 3*eps 602 za2=0.088_wp ! za2 = gam 603 za3=0.013_wp ! za3 = eps 604 ENDIF 605 606 zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) & 607 & + za2 * zsshb_e(:,:) + za3 * sshbb_e(:,:) 608 609 ! 610 ! Compute associated depths at U and V points: 611 IF ( lk_vvl.AND.(.NOT.ln_dynadv_vec) ) THEN 612 ! 613 DO jj = 2, jpjm1 614 DO ji = 2, jpim1 615 zx1 = z1_2 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) & 616 & * ( e1t(ji ,jj) * e2t(ji ,jj) * zsshp2_e(ji ,jj) & 617 & + e1t(ji+1,jj) * e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) 618 zy1 = z1_2 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) & 619 & * ( e1t(ji,jj ) * e2t(ji,jj ) * zsshp2_e(ji,jj ) & 620 & + e1t(ji,jj+1) * e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) 621 zhust_e(ji,jj) = hu_0(ji,jj) + zx1 622 zhvst_e(ji,jj) = hv_0(ji,jj) + zy1 623 END DO 624 END DO 625 ENDIF 626 ! 627 ! Add Coriolis trend: 628 ! zwz array below or triads normally depend on sea level with key_vvl and should be updated 629 ! at each time step. We however keep them constant here for optimization. 630 ! Recall that zwx and zwy arrays hold fluxes at this stage: 631 ! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5 632 ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 589 633 ! 590 634 IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN !== energy conserving or mixed scheme ==! … … 761 805 ! !* Time filter and swap 762 806 ! ! -------------------- 807 sshbb_e(:,:) = zsshb_e(:,:) 808 ubb_e (:,:) = zub_e (:,:) 809 vbb_e (:,:) = zvb_e (:,:) 763 810 IF( jn == 1 ) THEN ! Swap only (1st Euler time step) 764 811 zsshb_e(:,:) = sshn_e(:,:) … … 826 873 ! 827 874 ! !* Time average ==> after barotropic u, v, ssh 828 zcoef = 1._wp / ( 2 * nn_baro+ 1 )875 zcoef = 1._wp / ( icycle + 1 ) 829 876 zu_sum(:,:) = zcoef * zu_sum (:,:) 830 877 zv_sum(:,:) = zcoef * zv_sum (:,:) … … 855 902 CALL wrk_dealloc( jpi, jpj, zhu_b, zhv_b ) 856 903 CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zsshp2_e ) 904 CALL wrk_dealloc( jpi, jpj, zhust_e, zhvst_e ) 857 905 ! 858 906 IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_ts') … … 1052 1100 sshn_b(:,:) = sshb(:,:) ! if not in restart set previous time mean to current baroclinic before value 1053 1101 ENDIF 1102 1103 IF( .NOT.ln_bt_av .AND. iom_varid( numror, 'sshbb_e', ldstop = .FALSE. ) > 0) THEN 1104 CALL iom_get( numror, jpdom_autoglo, 'sshbb_e' , sshbb_e(:,:) ) 1105 CALL iom_get( numror, jpdom_autoglo, 'ubb_e' , ubb_e(:,:) ) 1106 CALL iom_get( numror, jpdom_autoglo, 'vbb_e' , vbb_e(:,:) ) 1107 CALL iom_get( numror, jpdom_autoglo, 'sshb_e' , sshb_e(:,:) ) 1108 CALL iom_get( numror, jpdom_autoglo, 'ub_e' , ub_e(:,:) ) 1109 CALL iom_get( numror, jpdom_autoglo, 'vb_e' , vb_e(:,:) ) 1110 ELSE 1111 sshbb_e = sshn_b ! ACC GUESS WORK 1112 ubb_e = ub_b 1113 vbb_e = vb_b 1114 sshb_e = sshn_b 1115 ub_e = ub_b 1116 vb_e = vb_b 1117 ENDIF 1054 1118 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 1055 1119 CALL iom_rstput( kt, nitrst, numrow, 'un_b' , un_b (:,:) ) ! external velocity and ssh 1056 1120 CALL iom_rstput( kt, nitrst, numrow, 'vn_b' , vn_b (:,:) ) ! issued from barotropic loop 1057 1121 CALL iom_rstput( kt, nitrst, numrow, 'sshn_b' , sshn_b(:,:) ) ! 1122 IF (.NOT.ln_bt_av) THEN 1123 CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:) ) 1124 CALL iom_rstput( kt, nitrst, numrow, 'ubb_e' , ubb_e(:,:) ) 1125 CALL iom_rstput( kt, nitrst, numrow, 'vbb_e' , vbb_e(:,:) ) 1126 CALL iom_rstput( kt, nitrst, numrow, 'sshb_e' , sshb_e(:,:) ) 1127 CALL iom_rstput( kt, nitrst, numrow, 'ub_e' , ub_e(:,:) ) 1128 CALL iom_rstput( kt, nitrst, numrow, 'vb_e' , vb_e(:,:) ) 1129 ENDIF 1058 1130 ENDIF 1059 1131 ! … … 1154 1226 1155 1227 #else 1156 !!---------------------------------------------------------------------- 1157 !! Default case : Empty module No standar t free surface cst volume1158 !!---------------------------------------------------------------------- 1228 !!--------------------------------------------------------------------------- 1229 !! Default case : Empty module No standard free surface constant volume 1230 !!--------------------------------------------------------------------------- 1159 1231 1160 1232 USE par_kind -
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
r4105 r4194 326 326 END DO 327 327 ! 328 DO jj = 2, jpjm1 !== th rid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk ==328 DO jj = 2, jpjm1 !== third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk == 329 329 DO ji = fs_2, fs_jpim1 ! vector opt. 330 330 va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)
Note: See TracChangeset
for help on using the changeset viewer.
