- Timestamp:
- 2015-12-16T10:25:22+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r5930 r6060 1 1 MODULE dynspg_ts 2 !!====================================================================== 3 !! *** MODULE dynspg_ts *** 4 !! Ocean dynamics: surface pressure gradient trend, split-explicit scheme 2 5 !!====================================================================== 3 6 !! History : 1.0 ! 2004-12 (L. Bessieres, G. Madec) Original code … … 13 16 !! 3.7 ! 2015-11 (J. Chanut) free surface simplification 14 17 !!--------------------------------------------------------------------- 18 15 19 !!---------------------------------------------------------------------- 16 !! split explicit free surface17 !! ----------------------------------------------------------------------18 !! dyn_spg_ts : compute surface pressure gradient trend using a time-19 !! splitting scheme and add to the general trend20 !! dyn_spg_ts : compute surface pressure gradient trend using a time-splitting scheme 21 !! dyn_spg_ts_init: initialisation of the time-splitting scheme 22 !! ts_wgt : set time-splitting weights for temporal averaging (or not) 23 !! ts_rst : read/write time-splitting fields in restart file 20 24 !!---------------------------------------------------------------------- 21 25 USE oce ! ocean dynamics and tracers 22 26 USE dom_oce ! ocean space and time domain 23 27 USE sbc_oce ! surface boundary condition: ocean 28 USE zdf_oce ! Bottom friction coefts 24 29 USE sbcisf ! ice shelf variable (fwfisf) 30 USE sbcapr ! surface boundary condition: atmospheric pressure 31 USE dynadv , ONLY: ln_dynadv_vec 25 32 USE phycst ! physical constants 26 33 USE dynvor ! vorticity term … … 30 37 USE sbctide ! tides 31 38 USE updtide ! tide potential 39 ! 40 USE in_out_manager ! I/O manager 32 41 USE lib_mpp ! distributed memory computing library 33 42 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 34 43 USE prtctl ! Print control 35 USE in_out_manager ! I/O manager36 44 USE iom ! IOM library 37 45 USE restart ! only for lrst_oce 38 USE zdf_oce ! Bottom friction coefts39 46 USE wrk_nemo ! Memory Allocation 40 47 USE timing ! Timing 41 USE sbcapr ! surface boundary condition: atmospheric pressure42 USE dynadv, ONLY: ln_dynadv_vec43 48 #if defined key_agrif 44 49 USE agrif_opa_interp ! agrif … … 59 64 REAL(wp),SAVE :: rdtbt ! Barotropic time step 60 65 61 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: & 62 wgtbtp1, & ! Primary weights used for time filtering of barotropic variables 63 wgtbtp2 ! Secondary weights used for time filtering of barotropic variables 64 65 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz ! ff/h at F points 66 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne ! triad of coriolis parameter 67 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse ! (only used with een vorticity scheme) 66 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: wgtbtp1, wgtbtp2 !: 1st & 2nd weights used in time filtering of barotropic fields 67 68 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz !: ff/h at F points 69 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne !: triad of coriolis parameter 70 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse !: (only used with een vorticity scheme) 71 72 !! Time filtered arrays at baroclinic time step: 73 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv , vn_adv !: Advection vel. at "now" barocl. step 68 74 69 75 !! * Substitutions 70 # include "domzgr_substitute.h90"71 76 # include "vectopt_loop_substitute.h90" 72 77 !!---------------------------------------------------------------------- … … 81 86 !! *** routine dyn_spg_ts_alloc *** 82 87 !!---------------------------------------------------------------------- 83 INTEGER :: ierr( 4)88 INTEGER :: ierr(3) 84 89 !!---------------------------------------------------------------------- 85 90 ierr(:) = 0 86 87 ALLOCATE( ssha_e(jpi,jpj), sshn_e(jpi,jpj), sshb_e(jpi,jpj), sshbb_e(jpi,jpj), & 88 & ua_e(jpi,jpj), un_e(jpi,jpj), ub_e(jpi,jpj), ubb_e(jpi,jpj), & 89 & va_e(jpi,jpj), vn_e(jpi,jpj), vb_e(jpi,jpj), vbb_e(jpi,jpj), & 90 & hu_e(jpi,jpj), hur_e(jpi,jpj), hv_e(jpi,jpj), hvr_e(jpi,jpj), STAT= ierr(1) ) 91 92 ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT= ierr(2) ) 93 94 IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 95 & ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) ) 96 97 ALLOCATE( ub2_b(jpi,jpj), vb2_b(jpi,jpj), un_adv(jpi,jpj), vn_adv(jpi,jpj), & 98 #if defined key_agrif 99 & ub2_i_b(jpi,jpj), vb2_i_b(jpi,jpj) , & 100 #endif 101 & STAT= ierr(4)) 102 103 dyn_spg_ts_alloc = MAXVAL(ierr(:)) 104 91 ! 92 ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT=ierr(1) ) 93 ! 94 IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 95 & ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(2) ) 96 ! 97 ALLOCATE( un_adv(jpi,jpj), vn_adv(jpi,jpj) , STAT=ierr(3) ) 98 ! 99 dyn_spg_ts_alloc = MAXVAL( ierr(:) ) 100 ! 105 101 IF( lk_mpp ) CALL mpp_sum( dyn_spg_ts_alloc ) 106 102 IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_warn('dyn_spg_ts_alloc: failed to allocate arrays') … … 112 108 !!---------------------------------------------------------------------- 113 109 !! 114 !! ** Purpose : 115 !! -Compute the now trend due to the explicit time stepping116 !! of the quasi-linear barotropic system.110 !! ** Purpose : - Compute the now trend due to the explicit time stepping 111 !! of the quasi-linear barotropic system, and add it to the 112 !! general momentum trend. 117 113 !! 118 !! ** Method : 114 !! ** Method : - split-explicit schem (time splitting) : 119 115 !! Barotropic variables are advanced from internal time steps 120 116 !! "n" to "n+1" if ln_bt_fw=T … … 130 126 !! continuity equation taken at the baroclinic time steps. This 131 127 !! ensures tracers conservation. 132 !! - Update 3d trend (ua, va)with barotropic component.128 !! - (ua, va) momentum trend updated with barotropic component. 133 129 !! 134 !! References : Shchepetkin, A.F. and J.C. McWilliams, 2005: 135 !! The regional oceanic modeling system (ROMS): 136 !! a split-explicit, free-surface, 137 !! topography-following-coordinate oceanic model. 138 !! Ocean Modelling, 9, 347-404. 130 !! References : Shchepetkin and McWilliams, Ocean Modelling, 2005. 139 131 !!--------------------------------------------------------------------- 140 !141 132 INTEGER, INTENT(in) :: kt ! ocean time-step index 142 133 ! 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 148 REAL(wp) :: zx1, zy1, zx2, zy2 149 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 ! - -150 REAL(wp) :: zu_spg, zv_spg 151 REAL(wp) :: zhura, zhvra ! - -152 REAL(wp) :: za0, za1, za2, za3 ! - -134 LOGICAL :: ll_fw_start ! if true, forward integration 135 LOGICAL :: ll_init ! if true, special startup of 2d equations 136 INTEGER :: ji, jj, jk, jn ! dummy loop indices 137 INTEGER :: ikbu, ikbv, noffset ! local integers 138 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars 139 REAL(wp) :: zx1, zy1, zx2, zy2 ! - - 140 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 ! - - 141 REAL(wp) :: zu_spg, zv_spg ! - - 142 REAL(wp) :: zhura, zhvra ! - - 143 REAL(wp) :: za0, za1, za2, za3 ! - - 153 144 ! 154 145 REAL(wp), POINTER, DIMENSION(:,:) :: zsshp2_e … … 163 154 ! 164 155 ! !* Allocate temporary arrays 165 CALL wrk_alloc( jpi, jpj, zsshp2_e, zhdiv ) 166 CALL wrk_alloc( jpi, jpj, zu_trd, zv_trd) 167 CALL wrk_alloc( jpi, jpj, zwx, zwy, 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 156 CALL wrk_alloc( jpi,jpj, zsshp2_e, zhdiv ) 157 CALL wrk_alloc( jpi,jpj, zu_trd, zv_trd) 158 CALL wrk_alloc( jpi,jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc) 159 CALL wrk_alloc( jpi,jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e) 160 CALL wrk_alloc( jpi,jpj, zsshu_a, zsshv_a ) 161 CALL wrk_alloc( jpi,jpj, zhf ) 162 ! 163 z1_12 = 1._wp / 12._wp !* Local constant initialization 174 164 z1_8 = 0.125_wp 175 165 z1_4 = 0.25_wp 176 166 z1_2 = 0.5_wp 177 167 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 168 ! ! reciprocal of baroclinic time step 169 IF( kt == nit000 .AND. neuler == 0 ) THEN ; z2dt_bf = rdt 170 ELSE ; z2dt_bf = 2.0_wp * rdt 183 171 ENDIF 184 172 z1_2dt_b = 1.0_wp / z2dt_bf 185 173 ! 186 ll_init = ln_bt_av! if no time averaging, then no specific restart174 ll_init = ln_bt_av ! if no time averaging, then no specific restart 187 175 ll_fw_start = .FALSE. 188 ! 189 ! time offset in steps for bdy data update 190 IF (.NOT.ln_bt_fw) THEN ; noffset=-nn_baro ; ELSE ; noffset = 0 ; ENDIF 176 ! ! time offset in steps for bdy data update 177 IF( .NOT.ln_bt_fw ) THEN ; noffset = - nn_baro 178 ELSE ; noffset = 0 179 ENDIF 191 180 ! 192 181 IF( kt == nit000 ) THEN !* initialisation … … 197 186 IF(lwp) WRITE(numout,*) 198 187 ! 199 IF (neuler==0)ll_init=.TRUE.200 ! 201 IF (ln_bt_fw.OR.(neuler==0)) THEN202 ll_fw_start=.TRUE.203 noffset= 0188 IF( neuler == 0 ) ll_init=.TRUE. 189 ! 190 IF( ln_bt_fw .OR. neuler == 0 ) THEN 191 ll_fw_start =.TRUE. 192 noffset = 0 204 193 ELSE 205 ll_fw_start=.FALSE.194 ll_fw_start =.FALSE. 206 195 ENDIF 207 196 ! 208 197 ! Set averaging weights and cycle length: 209 CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2) 210 ! 198 CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) 211 199 ! 212 200 ENDIF … … 219 207 ! and update depths at T-F points (ht and zhf resp.) at each barotropic time step 220 208 ! 221 IF ( kt == nit000 .OR. lk_vvl) THEN222 IF ( ln_dynvor_een ) THEN!== EEN scheme ==!209 IF( kt == nit000 .OR. .NOT.ln_linssh ) THEN 210 IF( ln_dynvor_een ) THEN !== EEN scheme ==! 223 211 SELECT CASE( nn_een_e3f ) !* ff/e3 at F-point 224 212 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 225 213 DO jj = 1, jpjm1 226 214 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._wp215 zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + & 216 & ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) * 0.25_wp 229 217 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff(ji,jj) / zwz(ji,jj) 230 218 END DO … … 233 221 DO jj = 1, jpjm1 234 222 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 ) ) &223 zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + & 224 & ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) & 237 225 & / ( MAX( 1._wp, tmask(ji ,jj+1, 1) + tmask(ji+1,jj+1, 1) + & 238 226 & tmask(ji ,jj , 1) + tmask(ji+1,jj , 1) ) ) … … 276 264 DO jk = 1, jpkm1 277 265 DO jj = 1, jpjm1 278 zhf(:,jj) = zhf(:,jj) + fse3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)266 zhf(:,jj) = zhf(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 279 267 END DO 280 268 END DO … … 308 296 ! 309 297 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)298 zu_frc(:,:) = zu_frc(:,:) + e3u_n(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 299 zv_frc(:,:) = zv_frc(:,:) + e3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 312 300 END DO 313 301 ! 314 zu_frc(:,:) = zu_frc(:,:) * hur(:,:)315 zv_frc(:,:) = zv_frc(:,:) * hvr(:,:)302 zu_frc(:,:) = zu_frc(:,:) * r1_hu_n(:,:) 303 zv_frc(:,:) = zv_frc(:,:) * r1_hv_n(:,:) 316 304 ! 317 305 ! … … 327 315 ! !* barotropic Coriolis trends (vorticity scheme dependent) 328 316 ! ! -------------------------------------------------------- 329 zwx(:,:) = un_b(:,:) * hu (:,:) * e2u(:,:) ! now fluxes330 zwy(:,:) = vn_b(:,:) * hv (:,:) * e1v(:,:)317 zwx(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:) ! now fluxes 318 zwy(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:) 331 319 ! 332 320 IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN ! energy conserving or mixed scheme … … 373 361 ! !* Right-Hand-Side of the barotropic momentum equation 374 362 ! ! ---------------------------------------------------- 375 IF( lk_vvl ) THEN! Variable volume : remove surface pressure gradient363 IF( .NOT.ln_linssh ) THEN ! Variable volume : remove surface pressure gradient 376 364 DO jj = 2, jpjm1 377 365 DO ji = fs_2, fs_jpim1 ! vector opt. … … 411 399 ! 412 400 ! 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(:,:)401 zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * bfrua(:,:) * zwx(:,:) 402 zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * bfrva(:,:) * zwy(:,:) 415 403 ! 416 IF (ln_bt_fw) THEN ! Add wind forcing417 zu_frc(:,:) = zu_frc(:,:) + zraur * utau(:,:) * hur(:,:)418 zv_frc(:,:) = zv_frc(:,:) + zraur * vtau(:,:) * hvr(:,:)404 IF( ln_bt_fw ) THEN ! Add wind forcing 405 zu_frc(:,:) = zu_frc(:,:) + zraur * utau(:,:) * r1_hu_n(:,:) 406 zv_frc(:,:) = zv_frc(:,:) + zraur * vtau(:,:) * r1_hv_n(:,:) 419 407 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(:,:)408 zu_frc(:,:) = zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * r1_hu_n(:,:) 409 zv_frc(:,:) = zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * r1_hv_n(:,:) 422 410 ENDIF 423 411 ! … … 484 472 ! 485 473 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields 486 sshn_e(:,:) = sshn(:,:)487 un_e (:,:) = un_b(:,:)488 vn_e (:,:) = vn_b(:,:)489 ! 490 hu_e (:,:) = hu(:,:)491 hv_e (:,:) = hv(:,:)492 hur_e (:,:) = hur(:,:)493 hvr_e (:,:) = hvr(:,:)474 sshn_e(:,:) = sshn(:,:) 475 un_e (:,:) = un_b(:,:) 476 vn_e (:,:) = vn_b(:,:) 477 ! 478 hu_e (:,:) = hu_n(:,:) 479 hv_e (:,:) = hv_n(:,:) 480 hur_e (:,:) = r1_hu_n(:,:) 481 hvr_e (:,:) = r1_hv_n(:,:) 494 482 ELSE ! CENTRED integration: start from BEFORE fields 495 sshn_e(:,:) = sshb(:,:)496 un_e (:,:) = ub_b(:,:)497 vn_e (:,:) = vb_b(:,:)498 ! 499 hu_e (:,:) = hu_b(:,:)500 hv_e (:,:) = hv_b(:,:)501 hur_e (:,:) = hur_b(:,:)502 hvr_e (:,:) = hvr_b(:,:)483 sshn_e(:,:) = sshb(:,:) 484 un_e (:,:) = ub_b(:,:) 485 vn_e (:,:) = vb_b(:,:) 486 ! 487 hu_e (:,:) = hu_b(:,:) 488 hv_e (:,:) = hv_b(:,:) 489 hur_e (:,:) = r1_hu_b(:,:) 490 hvr_e (:,:) = r1_hv_b(:,:) 503 491 ENDIF 504 492 ! … … 518 506 ! Update only tidal forcing at open boundaries 519 507 #if defined key_tide 520 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, time_offset=noffset)508 IF( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 ) 509 IF( ln_tide_pot .AND. lk_tide ) CALL upd_tide ( kt, kit=jn, time_offset= noffset ) 522 510 #endif 523 511 ! … … 537 525 va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:) 538 526 539 IF( lk_vvl ) THEN!* Update ocean depth (variable volume case only)527 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 540 528 ! ! ------------------ 541 529 ! Extrapolate Sea Level at step jit+0.5: … … 557 545 zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:) 558 546 ELSE 559 zhup2_e (:,:) = hu (:,:)560 zhvp2_e (:,:) = hv (:,:)547 zhup2_e (:,:) = hu_n(:,:) 548 zhvp2_e (:,:) = hv_n(:,:) 561 549 ENDIF 562 550 ! !* after ssh … … 569 557 ! 570 558 #if defined key_agrif 571 ! Set fluxes during predictor step to ensure 572 ! volume conservation 573 IF( (.NOT.Agrif_Root()).AND.ln_bt_fw ) THEN 559 ! Set fluxes during predictor step to ensure volume conservation 560 IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN 574 561 IF((nbondi == -1).OR.(nbondi == 2)) THEN 575 562 DO jj=1,jpj … … 611 598 612 599 #if defined key_bdy 613 ! Duplicate sea level across open boundaries (this is only cosmetic if l k_vvl=.false.)614 IF (lk_bdy)CALL bdy_ssh( ssha_e )600 ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T) 601 IF( lk_bdy ) CALL bdy_ssh( ssha_e ) 615 602 #endif 616 603 #if defined key_agrif 617 IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )604 IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn ) 618 605 #endif 619 606 ! 620 607 ! Sea Surface Height at u-,v-points (vvl case only) 621 IF ( lk_vvl) THEN608 IF( .NOT.ln_linssh ) THEN 622 609 DO jj = 2, jpjm1 623 610 DO ji = 2, jpim1 ! NO Vector Opt. … … 651 638 za3=0.013_wp ! za3 = eps 652 639 ENDIF 653 640 ! 654 641 zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) & 655 642 & + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 656 657 643 ! 658 644 ! Compute associated depths at U and V points: 659 IF ( lk_vvl.AND.(.NOT.ln_dynadv_vec) ) THEN645 IF( .NOT.ln_linssh .AND. .NOT.ln_dynadv_vec ) THEN !* Vector form 660 646 ! 661 647 DO jj = 2, jpjm1 … … 674 660 ! 675 661 ! Add Coriolis trend: 676 ! zwz array below or triads normally depend on sea level with key_vvland should be updated662 ! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated 677 663 ! at each time step. We however keep them constant here for optimization. 678 664 ! Recall that zwx and zwy arrays hold fluxes at this stage: … … 680 666 ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 681 667 ! 682 IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN 668 IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN !== energy conserving or mixed scheme ==! 683 669 DO jj = 2, jpjm1 684 670 DO ji = fs_2, fs_jpim1 ! vector opt. … … 692 678 END DO 693 679 ! 694 ELSEIF ( ln_dynvor_ens ) THEN 680 ELSEIF ( ln_dynvor_ens ) THEN !== enstrophy conserving scheme ==! 695 681 DO jj = 2, jpjm1 696 682 DO ji = fs_2, fs_jpim1 ! vector opt. … … 704 690 END DO 705 691 ! 706 ELSEIF ( ln_dynvor_een ) THEN !== energy and enstrophy conserving scheme ==!692 ELSEIF ( ln_dynvor_een ) THEN !== energy and enstrophy conserving scheme ==! 707 693 DO jj = 2, jpjm1 708 694 DO ji = fs_2, fs_jpim1 ! vector opt. … … 748 734 ! 749 735 ! Set next velocities: 750 IF( ln_dynadv_vec .OR. (.NOT. lk_vvl) ) THEN !Vector form736 IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form 751 737 DO jj = 2, jpjm1 752 738 DO ji = fs_2, fs_jpim1 ! vector opt. 753 ua_e(ji,jj) = ( un_e(ji,jj) & 754 & + rdtbt * ( zwx(ji,jj) & 755 & + zu_trd(ji,jj) & 756 & + zu_frc(ji,jj) ) & 757 & ) * umask(ji,jj,1) 758 759 va_e(ji,jj) = ( vn_e(ji,jj) & 760 & + rdtbt * ( zwy(ji,jj) & 761 & + zv_trd(ji,jj) & 762 & + zv_frc(ji,jj) ) & 763 & ) * vmask(ji,jj,1) 764 END DO 765 END DO 766 767 ELSE ! Flux form 739 ua_e(ji,jj) = ( un_e(ji,jj) & 740 & + rdtbt * ( zwx(ji,jj) & 741 & + zu_trd(ji,jj) & 742 & + zu_frc(ji,jj) ) ) * umask(ji,jj,1) 743 ! 744 va_e(ji,jj) = ( vn_e(ji,jj) & 745 & + rdtbt * ( zwy(ji,jj) & 746 & + zv_trd(ji,jj) & 747 & + zv_frc(ji,jj) ) ) * vmask(ji,jj,1) 748 END DO 749 END DO 750 ! 751 ELSE !* Flux form 768 752 DO jj = 2, jpjm1 769 753 DO ji = fs_2, fs_jpim1 ! vector opt. 770 771 zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1)) 772 zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1)) 773 774 ua_e(ji,jj) = ( hu_e(ji,jj) * un_e(ji,jj) & 775 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 776 & + zhup2_e(ji,jj) * zu_trd(ji,jj) & 777 & + hu(ji,jj) * zu_frc(ji,jj) ) & 778 & ) * zhura 779 780 va_e(ji,jj) = ( hv_e(ji,jj) * vn_e(ji,jj) & 781 & + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) & 782 & + zhvp2_e(ji,jj) * zv_trd(ji,jj) & 783 & + hv(ji,jj) * zv_frc(ji,jj) ) & 784 & ) * zhvra 785 END DO 786 END DO 787 ENDIF 788 ! 789 IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only) 790 ! ! ---------------------------------------------- 754 zhura = umask(ji,jj,1) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1) ) 755 zhvra = vmask(ji,jj,1) / ( hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1) ) 756 ! 757 ua_e(ji,jj) = ( hu_e(ji,jj) * un_e(ji,jj) & 758 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 759 & + zhup2_e(ji,jj) * zu_trd(ji,jj) & 760 & + hu_n(ji,jj) * zu_frc(ji,jj) ) ) * zhura 761 ! 762 va_e(ji,jj) = ( hv_e(ji,jj) * vn_e(ji,jj) & 763 & + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) & 764 & + zhvp2_e(ji,jj) * zv_trd(ji,jj) & 765 & + hv_n(ji,jj) * zv_frc(ji,jj) ) ) * zhvra 766 END DO 767 END DO 768 ENDIF 769 ! 770 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 791 771 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 792 772 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) … … 795 775 ! 796 776 ENDIF 797 ! !* domain lateral boundary 798 ! ! ----------------------- 799 ! 777 ! !* domain lateral boundary 800 778 CALL lbc_lnk_multi( ua_e, 'U', -1._wp, va_e , 'V', -1._wp ) 801 779 ! 802 780 #if defined key_bdy 803 804 IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e )781 ! ! open boundaries 782 IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e ) 805 783 #endif 806 784 #if defined key_agrif … … 824 802 ! ! ---------------------- 825 803 za1 = wgtbtp1(jn) 826 IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN ! Sum velocities804 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! Sum velocities 827 805 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) 828 806 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) … … 843 821 zwx(:,:) = un_adv(:,:) 844 822 zwy(:,:) = vn_adv(:,:) 845 IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN846 un_adv(:,:) = zwx(:,:) *hur(:,:)847 vn_adv(:,:) = zwy(:,:) *hvr(:,:)823 IF( ( kt == nit000 .AND. neuler==0 ) .OR. .NOT.ln_bt_fw ) THEN 824 un_adv(:,:) = zwx(:,:) * r1_hu_n(:,:) 825 vn_adv(:,:) = zwy(:,:) * r1_hv_n(:,:) 848 826 ELSE 849 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:) ) * hur(:,:)850 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zwy(:,:) ) * hvr(:,:)827 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:) ) * r1_hu_n(:,:) 828 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zwy(:,:) ) * r1_hv_n(:,:) 851 829 END IF 852 830 853 IF (ln_bt_fw) THEN ! Save integrated transport for next computation831 IF( ln_bt_fw ) THEN ! Save integrated transport for next computation 854 832 ub2_b(:,:) = zwx(:,:) 855 833 vb2_b(:,:) = zwy(:,:) … … 857 835 ! 858 836 ! Update barotropic trend: 859 IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN837 IF( ln_dynadv_vec .OR. ln_linssh ) THEN 860 838 DO jk=1,jpkm1 861 839 ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b … … 877 855 ! 878 856 DO jk=1,jpkm1 879 ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b880 va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b857 ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b 858 va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b 881 859 END DO 882 860 ! Save barotropic velocities not transport: 883 ua_b 884 va_b 861 ua_b(:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - umask(:,:,1) ) 862 va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - vmask(:,:,1) ) 885 863 ENDIF 886 864 ! 887 865 DO jk = 1, jpkm1 888 866 ! Correct velocities: 889 un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) ) *umask(:,:,jk)890 vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) ) *vmask(:,:,jk)867 un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) ) * umask(:,:,jk) 868 vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) ) * vmask(:,:,jk) 891 869 ! 892 870 END DO 871 ! 872 CALL iom_put( "ubar", un_adv(:,:) ) ! barotropic i-current 873 CALL iom_put( "vbar", vn_adv(:,:) ) ! barotropic i-current 893 874 ! 894 875 #if defined key_agrif … … 896 877 ! (used to update coarse grid transports at next time step) 897 878 ! 898 IF ( (.NOT.Agrif_Root()).AND.(ln_bt_fw)) THEN899 IF ( Agrif_NbStepint().EQ.0 ) THEN900 ub2_i_b(:,:) = 0. e0901 vb2_i_b(:,:) = 0. e0879 IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN 880 IF( Agrif_NbStepint() == 0 ) THEN 881 ub2_i_b(:,:) = 0._wp 882 vb2_i_b(:,:) = 0._wp 902 883 END IF 903 884 ! … … 906 887 vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:) 907 888 ENDIF 908 !909 !910 889 #endif 911 !912 890 ! !* write time-spliting arrays in the restart 913 IF( lrst_oce .AND.ln_bt_fw) CALL ts_rst( kt, 'WRITE' )914 ! 915 CALL wrk_dealloc( jpi, jpj,zsshp2_e, zhdiv )916 CALL wrk_dealloc( jpi, jpj,zu_trd, zv_trd )917 CALL wrk_dealloc( jpi, jpj,zwx, zwy, zssh_frc, zu_frc, zv_frc )918 CALL wrk_dealloc( jpi, jpj,zhup2_e, zhvp2_e, zhust_e, zhvst_e )919 CALL wrk_dealloc( jpi, jpj,zsshu_a, zsshv_a )920 CALL wrk_dealloc( jpi, jpj,zhf )891 IF( lrst_oce .AND.ln_bt_fw ) CALL ts_rst( kt, 'WRITE' ) 892 ! 893 CALL wrk_dealloc( jpi,jpj, zsshp2_e, zhdiv ) 894 CALL wrk_dealloc( jpi,jpj, zu_trd, zv_trd ) 895 CALL wrk_dealloc( jpi,jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc ) 896 CALL wrk_dealloc( jpi,jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e ) 897 CALL wrk_dealloc( jpi,jpj, zsshu_a, zsshv_a ) 898 CALL wrk_dealloc( jpi,jpj, zhf ) 921 899 ! 922 900 IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_ts') 923 901 ! 924 902 END SUBROUTINE dyn_spg_ts 903 925 904 926 905 SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2) … … 1001 980 END SUBROUTINE ts_wgt 1002 981 982 1003 983 SUBROUTINE ts_rst( kt, cdrw ) 1004 984 !!--------------------------------------------------------------------- … … 1054 1034 END SUBROUTINE ts_rst 1055 1035 1056 SUBROUTINE dyn_spg_ts_init( kt ) 1036 1037 SUBROUTINE dyn_spg_ts_init 1057 1038 !!--------------------------------------------------------------------- 1058 1039 !! *** ROUTINE dyn_spg_ts_init *** … … 1060 1041 !! ** Purpose : Set time splitting options 1061 1042 !!---------------------------------------------------------------------- 1062 INTEGER , INTENT(in) :: kt ! ocean time-step 1063 ! 1064 INTEGER :: ji ,jj 1065 REAL(wp) :: zxr2, zyr2, zcmax 1066 REAL(wp), POINTER, DIMENSION(:,:) :: zcu 1067 !! 1043 INTEGER :: ji ,jj ! dummy loop indices 1044 REAL(wp) :: zxr2, zyr2, zcmax ! local scalar 1045 REAL(wp), POINTER, DIMENSION(:,:) :: zcu 1068 1046 !!---------------------------------------------------------------------- 1069 1047 ! 1070 1048 ! Max courant number for ext. grav. waves 1071 1049 ! 1072 CALL wrk_alloc( jpi, jpj,zcu )1050 CALL wrk_alloc( jpi,jpj, zcu ) 1073 1051 ! 1074 1052 DO jj = 1, jpj … … 1084 1062 1085 1063 ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax 1086 IF (ln_bt_auto)nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)1064 IF( ln_bt_auto ) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax) 1087 1065 1088 1066 rdtbt = rdt / REAL( nn_baro , wp ) … … 1114 1092 #if defined key_agrif 1115 1093 ! Restrict the use of Agrif to the forward case only 1116 IF ((.NOT.ln_bt_fw ).AND.(.NOT.Agrif_Root()))CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )1094 IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' ) 1117 1095 #endif 1118 1096 ! 1119 1097 IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt 1120 1098 SELECT CASE ( nn_bt_flt ) 1121 CASE( 0 ); IF(lwp) WRITE(numout,*) ' Dirac'1122 CASE( 1 ); IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_baro'1123 CASE( 2 ); IF(lwp) WRITE(numout,*) ' Boxcar: width = 2*nn_baro'1124 1099 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac' 1100 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_baro' 1101 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = 2*nn_baro' 1102 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' ) 1125 1103 END SELECT 1126 1104 ! … … 1130 1108 IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax 1131 1109 ! 1132 IF ((.NOT.ln_bt_av).AND.(.NOT.ln_bt_fw)) THEN1110 IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN 1133 1111 CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' ) 1134 1112 ENDIF 1135 IF 1113 IF( zcmax>0.9_wp ) THEN 1136 1114 CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' ) 1137 1115 ENDIF 1138 1116 ! 1139 CALL wrk_dealloc( jpi, jpj,zcu )1117 CALL wrk_dealloc( jpi,jpj, zcu ) 1140 1118 ! 1141 1119 END SUBROUTINE dyn_spg_ts_init
Note: See TracChangeset
for help on using the changeset viewer.
