Changeset 5845 for branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
- Timestamp:
- 2015-10-31T08:40:45+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r5836 r5845 75 75 76 76 !! * Substitutions 77 # include "domzgr_substitute.h90"78 77 # include "vectopt_loop_substitute.h90" 79 78 !!---------------------------------------------------------------------- … … 91 90 !!---------------------------------------------------------------------- 92 91 ierr(:) = 0 93 92 ! 94 93 ALLOCATE( sshb_e(jpi,jpj), sshbb_e(jpi,jpj), & 95 94 & ub_e(jpi,jpj) , vb_e(jpi,jpj) , & 96 95 & ubb_e(jpi,jpj) , vbb_e(jpi,jpj) , STAT= ierr(1) ) 97 96 ! 98 97 ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT= ierr(2) ) 99 98 ! 100 99 IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 101 100 & ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) ) 102 103 dyn_spg_ts_alloc = MAXVAL(ierr(:)) 104 101 ! 102 dyn_spg_ts_alloc = MAXVAL( ierr(:) ) 105 103 IF( lk_mpp ) CALL mpp_sum( dyn_spg_ts_alloc ) 106 104 IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_warn('dynspg_oce_alloc: failed to allocate arrays') … … 138 136 !! Ocean Modelling, 9, 347-404. 139 137 !!--------------------------------------------------------------------- 140 !141 138 INTEGER, INTENT(in) :: kt ! ocean time-step index 142 139 ! 143 LOGICAL :: ll_fw_start 144 LOGICAL :: ll_init 145 INTEGER :: ji, jj, jk, jn 146 INTEGER :: ikbu, ikbv, noffset ! local integers147 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf 140 LOGICAL :: ll_fw_start ! if true, forward integration 141 LOGICAL :: ll_init ! if true, special startup of 2d equations 142 INTEGER :: ji, jj, jk, jn ! dummy loop indices 143 INTEGER :: ikbu, ikbv, noffset ! local integers 144 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars 148 145 REAL(wp) :: zx1, zy1, zx2, zy2 ! - - 149 146 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 ! - - 150 147 REAL(wp) :: zu_spg, zv_spg ! - - 151 148 REAL(wp) :: zhura, zhvra ! - - 152 REAL(wp) :: za0, za1, za2, za3 153 ! 154 REAL(wp), POINTER, DIMENSION(:,:) :: zun_e, zvn_e, zsshp2_e155 REAL(wp), POINTER, DIMENSION(:,:) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc156 REAL(wp), POINTER, DIMENSION(:,:) :: zu_sum, zv_sum, zwx, zwy, zhdiv157 REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e158 REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a159 REAL(wp), POINTER, DIMENSION(:,:) :: zhf149 REAL(wp) :: za0, za1, za2, za3 ! - - 150 ! 151 REAL(wp), POINTER, DIMENSION(:,:) :: zun_e, zvn_e, zsshp2_e 152 REAL(wp), POINTER, DIMENSION(:,:) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc 153 REAL(wp), POINTER, DIMENSION(:,:) :: zu_sum, zv_sum, zwx, zwy, zhdiv 154 REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e 155 REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a 156 REAL(wp), POINTER, DIMENSION(:,:) :: zhf 160 157 !!---------------------------------------------------------------------- 161 158 ! 162 159 IF( nn_timing == 1 ) CALL timing_start('dyn_spg_ts') 163 160 ! 164 ! !* Allocate temporary arrays 165 CALL wrk_alloc( jpi, jpj, zsshp2_e, zhdiv ) 166 CALL wrk_alloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e ) 167 CALL wrk_alloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc) 168 CALL wrk_alloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e) 169 CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a ) 170 CALL wrk_alloc( jpi, jpj, zhf ) 171 ! 172 ! !* Local constant initialization 173 z1_12 = 1._wp / 12._wp 161 CALL wrk_alloc( jpi,jpj, zsshp2_e, zhdiv ) 162 CALL wrk_alloc( jpi,jpj, zu_trd, zv_trd, zun_e, zvn_e ) 163 CALL wrk_alloc( jpi,jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc) 164 CALL wrk_alloc( jpi,jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e) 165 CALL wrk_alloc( jpi,jpj, zsshu_a, zsshv_a ) 166 CALL wrk_alloc( jpi,jpj, zhf ) 167 ! 168 z1_12 = 1._wp / 12._wp !* Local constant initialization 174 169 z1_8 = 0.125_wp 175 170 z1_4 = 0.25_wp 176 171 z1_2 = 0.5_wp 177 172 zraur = 1._wp / rau0 178 ! 179 IF( kt == nit000 .AND. neuler == 0 ) THEN ! reciprocal of baroclinic time step 180 z2dt_bf = rdt 181 ELSE 182 z2dt_bf = 2.0_wp * rdt 173 ! ! reciprocal of baroclinic time step 174 IF( kt == nit000 .AND. neuler == 0 ) THEN ; z2dt_bf = rdt 175 ELSE ; z2dt_bf = 2.0_wp * rdt 183 176 ENDIF 184 177 z1_2dt_b = 1.0_wp / z2dt_bf 185 178 ! 186 ll_init = ln_bt_av! if no time averaging, then no specific restart179 ll_init = ln_bt_av ! if no time averaging, then no specific restart 187 180 ll_fw_start = .FALSE. 188 ! 189 ! time offset in steps for bdy data update 190 IF (.NOT.ln_bt_fw) THEN ; noffset=-2*nn_baro ; ELSE ; noffset = 0 ; ENDIF 181 ! ! time offset in steps for bdy data update 182 IF( .NOT.ln_bt_fw ) THEN ; noffset =-2*nn_baro 183 ELSE ; noffset = 0 184 ENDIF 191 185 ! 192 186 IF( kt == nit000 ) THEN !* initialisation … … 197 191 IF(lwp) WRITE(numout,*) 198 192 ! 199 IF (neuler==0)ll_init=.TRUE.200 ! 201 IF (ln_bt_fw.OR.(neuler==0)) THEN202 ll_fw_start=.TRUE.203 noffset = 0193 IF( neuler == 0 ) ll_init=.TRUE. 194 ! 195 IF( ln_bt_fw .OR. neuler == 0 ) THEN 196 ll_fw_start=.TRUE. 197 noffset = 0 204 198 ELSE 205 ll_fw_start=.FALSE.199 ll_fw_start=.FALSE. 206 200 ENDIF 207 201 ! 208 202 ! Set averaging weights and cycle length: 209 CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2) 210 ! 203 CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) 211 204 ! 212 205 ENDIF … … 225 218 DO jj = 1, jpjm1 226 219 DO ji = 1, jpim1 227 zwz(ji,jj) = ( ht (ji ,jj+1) + ht(ji+1,jj+1) + &228 & ht (ji ,jj ) + ht(ji+1,jj ) ) / 4._wp220 zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + & 221 & ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) / 4._wp 229 222 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff(ji,jj) / zwz(ji,jj) 230 223 END DO … … 233 226 DO jj = 1, jpjm1 234 227 DO ji = 1, jpim1 235 zwz(ji,jj) = ( ht (ji ,jj+1) + ht(ji+1,jj+1) + &236 & ht (ji ,jj ) + ht(ji+1,jj ) ) &228 zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + & 229 & ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) & 237 230 & / ( MAX( 1._wp, tmask(ji ,jj+1, 1) + tmask(ji+1,jj+1, 1) + & 238 231 & tmask(ji ,jj , 1) + tmask(ji+1,jj , 1) ) ) … … 276 269 DO jk = 1, jpkm1 277 270 DO jj = 1, jpjm1 278 zhf(:,jj) = zhf(:,jj) + fse3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)271 zhf(:,jj) = zhf(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 279 272 END DO 280 273 END DO … … 308 301 ! 309 302 DO jk = 1, jpkm1 310 zu_frc(:,:) = zu_frc(:,:) + fse3u_n(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)311 zv_frc(:,:) = zv_frc(:,:) + fse3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)303 zu_frc(:,:) = zu_frc(:,:) + e3u_n(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 304 zv_frc(:,:) = zv_frc(:,:) + e3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 312 305 END DO 313 306 ! 314 zu_frc(:,:) = zu_frc(:,:) * hur(:,:)315 zv_frc(:,:) = zv_frc(:,:) * hvr(:,:)307 zu_frc(:,:) = zu_frc(:,:) * r1_hu_n(:,:) 308 zv_frc(:,:) = zv_frc(:,:) * r1_hv_n(:,:) 316 309 ! 317 310 ! … … 327 320 ! !* barotropic Coriolis trends (vorticity scheme dependent) 328 321 ! ! -------------------------------------------------------- 329 zwx(:,:) = un_b(:,:) * hu (:,:) * e2u(:,:) ! now fluxes330 zwy(:,:) = vn_b(:,:) * hv (:,:) * e1v(:,:)322 zwx(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:) ! now fluxes 323 zwy(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:) 331 324 ! 332 325 IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN ! energy conserving or mixed scheme … … 411 404 ! 412 405 ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 413 zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:)414 zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:)406 zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * bfrua(:,:) * zwx(:,:) 407 zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * bfrva(:,:) * zwy(:,:) 415 408 ! 416 409 IF (ln_bt_fw) THEN ! Add wind forcing 417 zu_frc(:,:) = zu_frc(:,:) + zraur * utau(:,:) * hur(:,:)418 zv_frc(:,:) = zv_frc(:,:) + zraur * vtau(:,:) * hvr(:,:)410 zu_frc(:,:) = zu_frc(:,:) + zraur * utau(:,:) * r1_hu_n(:,:) 411 zv_frc(:,:) = zv_frc(:,:) + zraur * vtau(:,:) * r1_hv_n(:,:) 419 412 ELSE 420 zu_frc(:,:) = zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * hur(:,:)421 zv_frc(:,:) = zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * hvr(:,:)413 zu_frc(:,:) = zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * r1_hu_n(:,:) 414 zv_frc(:,:) = zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * r1_hv_n(:,:) 422 415 ENDIF 423 416 ! … … 484 477 ! 485 478 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields 486 sshn_e(:,:) = sshn 487 zun_e (:,:) = un_b 488 zvn_e (:,:) = vn_b 489 ! 490 hu_e (:,:) = hu(:,:)491 hv_e (:,:) = hv(:,:)492 hur_e (:,:) = hur(:,:)493 hvr_e (:,:) = hvr(:,:)479 sshn_e(:,:) = sshn(:,:) 480 zun_e (:,:) = un_b(:,:) 481 zvn_e (:,:) = vn_b(:,:) 482 ! 483 hu_e (:,:) = hu_n(:,:) 484 hv_e (:,:) = hv_n(:,:) 485 hur_e (:,:) = r1_hu_n(:,:) 486 hvr_e (:,:) = r1_hv_n(:,:) 494 487 ELSE ! CENTRED integration: start from BEFORE fields 495 sshn_e(:,:) = sshb 496 zun_e (:,:) = ub_b 497 zvn_e (:,:) = vb_b 498 ! 499 hu_e (:,:) = hu_b(:,:)500 hv_e (:,:) = hv_b(:,:)501 hur_e (:,:) = hur_b(:,:)502 hvr_e (:,:) = hvr_b(:,:)488 sshn_e(:,:) = sshb(:,:) 489 zun_e (:,:) = ub_b(:,:) 490 zvn_e (:,:) = vb_b(:,:) 491 ! 492 hu_e (:,:) = hu_b(:,:) 493 hv_e (:,:) = hv_b(:,:) 494 hur_e (:,:) = r1_hu_b(:,:) 495 hvr_e (:,:) = r1_hv_b(:,:) 503 496 ENDIF 504 497 ! … … 519 512 #if defined key_tide 520 513 IF ( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1) ) 521 IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide ( kt, kit=jn, koffset=noffset )514 IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide ( kt, kit=jn, koffset=noffset ) 522 515 #endif 523 516 ! … … 557 550 zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:) 558 551 ELSE 559 zhup2_e (:,:) = hu (:,:)560 zhvp2_e (:,:) = hv (:,:)552 zhup2_e (:,:) = hu_n(:,:) 553 zhvp2_e (:,:) = hv_n(:,:) 561 554 ENDIF 562 555 ! !* after ssh … … 775 768 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 776 769 & + zhup2_e(ji,jj) * zu_trd(ji,jj) & 777 & + hu(ji,jj) * zu_frc(ji,jj) ) &770 & + hu_n(ji,jj) * zu_frc(ji,jj) ) & 778 771 & ) * zhura 779 772 … … 781 774 & + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) & 782 775 & + zhvp2_e(ji,jj) * zv_trd(ji,jj) & 783 & + hv(ji,jj) * zv_frc(ji,jj) ) &776 & + hv_n(ji,jj) * zv_frc(ji,jj) ) & 784 777 & ) * zhvra 785 778 END DO … … 857 850 ! 858 851 ! Set advection velocity correction: 859 IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN860 un_adv(:,:) = zu_sum(:,:) *hur(:,:)861 vn_adv(:,:) = zv_sum(:,:) *hvr(:,:)852 IF( ( kt == nit000 .AND. neuler==0 ) .OR. .NOT.ln_bt_fw ) THEN 853 un_adv(:,:) = zu_sum(:,:) * r1_hu_n(:,:) 854 vn_adv(:,:) = zv_sum(:,:) * r1_hv_n(:,:) 862 855 ELSE 863 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zu_sum(:,:) ) * hur(:,:)864 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zv_sum(:,:) ) * hvr(:,:)856 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zu_sum(:,:) ) * r1_hu_n(:,:) 857 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zv_sum(:,:) ) * r1_hv_n(:,:) 865 858 END IF 866 859 867 IF (ln_bt_fw) THEN ! Save integrated transport for next computation860 IF( ln_bt_fw ) THEN ! Save integrated transport for next computation 868 861 ub2_b(:,:) = zu_sum(:,:) 869 862 vb2_b(:,:) = zv_sum(:,:) … … 871 864 ! 872 865 ! Update barotropic trend: 873 IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN866 IF( ln_dynadv_vec .OR. .NOT.lk_vvl ) THEN 874 867 DO jk=1,jpkm1 875 868 ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b … … 878 871 ELSE 879 872 DO jk=1,jpkm1 880 ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b881 va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b873 ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b 874 va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b 882 875 END DO 883 876 ! Save barotropic velocities not transport: 884 ua_b 885 va_b 877 ua_b(:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - umask(:,:,1) ) 878 va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - vmask(:,:,1) ) 886 879 ENDIF 887 880 ! 888 881 DO jk = 1, jpkm1 889 882 ! Correct velocities: 890 un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) ) *umask(:,:,jk)891 vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) ) *vmask(:,:,jk)883 un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) ) * umask(:,:,jk) 884 vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) ) * vmask(:,:,jk) 892 885 ! 893 886 END DO … … 897 890 ! (used to update coarse grid transports at next time step) 898 891 ! 899 IF ( (.NOT.Agrif_Root()).AND.(ln_bt_fw)) THEN900 IF ( Agrif_NbStepint().EQ.0 ) THEN901 ub2_i_b(:,:) = 0. e0902 vb2_i_b(:,:) = 0. e0892 IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN 893 IF( Agrif_NbStepint() == 0 ) THEN 894 ub2_i_b(:,:) = 0._wp 895 vb2_i_b(:,:) = 0._wp 903 896 END IF 904 897 ! … … 912 905 ! 913 906 ! !* write time-spliting arrays in the restart 914 IF( lrst_oce .AND.ln_bt_fw) CALL ts_rst( kt, 'WRITE' )915 ! 916 CALL wrk_dealloc( jpi, jpj,zsshp2_e, zhdiv )917 CALL wrk_dealloc( jpi, jpj,zu_trd, zv_trd, zun_e, zvn_e )918 CALL wrk_dealloc( jpi, jpj,zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc )919 CALL wrk_dealloc( jpi, jpj,zhup2_e, zhvp2_e, zhust_e, zhvst_e )920 CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a)921 CALL wrk_dealloc( jpi, jpj,zhf )907 IF( lrst_oce .AND.ln_bt_fw ) CALL ts_rst( kt, 'WRITE' ) 908 ! 909 CALL wrk_dealloc( jpi,jpj, zsshp2_e, zhdiv ) 910 CALL wrk_dealloc( jpi,jpj, zu_trd, zv_trd, zun_e, zvn_e ) 911 CALL wrk_dealloc( jpi,jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc ) 912 CALL wrk_dealloc( jpi,jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e ) 913 CALL wrk_dealloc( jpi,jpj, zsshu_a, zsshv_a ) 914 CALL wrk_dealloc( jpi,jpj, zhf ) 922 915 ! 923 916 IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_ts') 924 917 ! 925 918 END SUBROUTINE dyn_spg_ts 919 926 920 927 921 SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2) … … 1094 1088 END DO 1095 1089 ELSE 1090 !!gm BUG ?? restartability issue if ssh changes are large.... 1091 !!gm We should just test this with ht_0 only, no? 1096 1092 DO jj = 1, jpj 1097 1093 DO ji =1, jpi 1098 1094 zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 1099 1095 zyr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 1100 zcu(ji,jj) = SQRT( grav * ht (ji,jj) * (zxr2 + zyr2) )1096 zcu(ji,jj) = SQRT( grav * ht_n(ji,jj) * (zxr2 + zyr2) ) 1101 1097 END DO 1102 1098 END DO
Note: See TracChangeset
for help on using the changeset viewer.
