- Timestamp:
- 2016-11-28T17:04:10+01:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r5930 r7351 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 34 USE wet_dry ! wetting/drying flux limter 27 35 USE bdy_par ! for lk_bdy 28 36 USE bdytides ! open boundary condition data … … 30 38 USE sbctide ! tides 31 39 USE updtide ! tide potential 40 ! 41 USE in_out_manager ! I/O manager 32 42 USE lib_mpp ! distributed memory computing library 33 43 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 34 44 USE prtctl ! Print control 35 USE in_out_manager ! I/O manager36 45 USE iom ! IOM library 37 46 USE restart ! only for lrst_oce 38 USE zdf_oce ! Bottom friction coefts39 47 USE wrk_nemo ! Memory Allocation 40 48 USE timing ! Timing 41 USE sbcapr ! surface boundary condition: atmospheric pressure 42 USE dynadv, ONLY: ln_dynadv_vec 49 USE diatmb ! Top,middle,bottom output 43 50 #if defined key_agrif 44 51 USE agrif_opa_interp ! agrif … … 48 55 #endif 49 56 57 50 58 IMPLICIT NONE 51 59 PRIVATE … … 59 67 REAL(wp),SAVE :: rdtbt ! Barotropic time step 60 68 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) 69 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: wgtbtp1, wgtbtp2 !: 1st & 2nd weights used in time filtering of barotropic fields 70 71 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz !: ff/h at F points 72 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne !: triad of coriolis parameter 73 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse !: (only used with een vorticity scheme) 74 75 !! Time filtered arrays at baroclinic time step: 76 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv , vn_adv !: Advection vel. at "now" barocl. step 68 77 69 78 !! * Substitutions 70 # include "domzgr_substitute.h90"71 79 # include "vectopt_loop_substitute.h90" 72 80 !!---------------------------------------------------------------------- … … 81 89 !! *** routine dyn_spg_ts_alloc *** 82 90 !!---------------------------------------------------------------------- 83 INTEGER :: ierr( 4)91 INTEGER :: ierr(3) 84 92 !!---------------------------------------------------------------------- 85 93 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 94 ! 95 ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT=ierr(1) ) 96 ! 97 IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 98 & ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(2) ) 99 ! 100 ALLOCATE( un_adv(jpi,jpj), vn_adv(jpi,jpj) , STAT=ierr(3) ) 101 ! 102 dyn_spg_ts_alloc = MAXVAL( ierr(:) ) 103 ! 105 104 IF( lk_mpp ) CALL mpp_sum( dyn_spg_ts_alloc ) 106 105 IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_warn('dyn_spg_ts_alloc: failed to allocate arrays') … … 112 111 !!---------------------------------------------------------------------- 113 112 !! 114 !! ** Purpose : 115 !! -Compute the now trend due to the explicit time stepping116 !! of the quasi-linear barotropic system.113 !! ** Purpose : - Compute the now trend due to the explicit time stepping 114 !! of the quasi-linear barotropic system, and add it to the 115 !! general momentum trend. 117 116 !! 118 !! ** Method : 117 !! ** Method : - split-explicit schem (time splitting) : 119 118 !! Barotropic variables are advanced from internal time steps 120 119 !! "n" to "n+1" if ln_bt_fw=T … … 130 129 !! continuity equation taken at the baroclinic time steps. This 131 130 !! ensures tracers conservation. 132 !! - Update 3d trend (ua, va)with barotropic component.131 !! - (ua, va) momentum trend updated with barotropic component. 133 132 !! 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. 133 !! References : Shchepetkin and McWilliams, Ocean Modelling, 2005. 139 134 !!--------------------------------------------------------------------- 140 !141 135 INTEGER, INTENT(in) :: kt ! ocean time-step index 142 136 ! 143 137 LOGICAL :: ll_fw_start ! if true, forward integration 144 138 LOGICAL :: ll_init ! if true, special startup of 2d equations 139 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables used in W/D 145 140 INTEGER :: ji, jj, jk, jn ! dummy loop indices 146 141 INTEGER :: ikbu, ikbv, noffset ! local integers 142 INTEGER :: iktu, iktv ! local integers 143 REAL(wp) :: zmdi 147 144 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars 148 145 REAL(wp) :: zx1, zy1, zx2, zy2 ! - - … … 158 155 REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a 159 156 REAL(wp), POINTER, DIMENSION(:,:) :: zhf 157 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef. 158 REAL(wp), POINTER, DIMENSION(:,:) :: wduflt1, wdvflt1 ! Wetting/Dying velocity filter coef. 160 159 !!---------------------------------------------------------------------- 161 160 ! … … 163 162 ! 164 163 ! !* 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 ! 164 CALL wrk_alloc( jpi,jpj, zsshp2_e, zhdiv ) 165 CALL wrk_alloc( jpi,jpj, zu_trd, zv_trd) 166 CALL wrk_alloc( jpi,jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc) 167 CALL wrk_alloc( jpi,jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e) 168 CALL wrk_alloc( jpi,jpj, zsshu_a, zsshv_a ) 169 CALL wrk_alloc( jpi,jpj, zhf ) 170 IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 ) 171 ! 172 zmdi=1.e+20 ! missing data indicator for masking 172 173 ! !* Local constant initialization 173 174 z1_12 = 1._wp / 12._wp … … 176 177 z1_2 = 0.5_wp 177 178 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 179 ! ! reciprocal of baroclinic time step 180 IF( kt == nit000 .AND. neuler == 0 ) THEN ; z2dt_bf = rdt 181 ELSE ; z2dt_bf = 2.0_wp * rdt 183 182 ENDIF 184 183 z1_2dt_b = 1.0_wp / z2dt_bf 185 184 ! 186 ll_init = ln_bt_av! if no time averaging, then no specific restart185 ll_init = ln_bt_av ! if no time averaging, then no specific restart 187 186 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 187 ! ! time offset in steps for bdy data update 188 IF( .NOT.ln_bt_fw ) THEN ; noffset = - nn_baro 189 ELSE ; noffset = 0 190 ENDIF 191 191 ! 192 192 IF( kt == nit000 ) THEN !* initialisation … … 197 197 IF(lwp) WRITE(numout,*) 198 198 ! 199 IF (neuler==0)ll_init=.TRUE.200 ! 201 IF (ln_bt_fw.OR.(neuler==0)) THEN202 ll_fw_start=.TRUE.203 noffset= 0199 IF( neuler == 0 ) ll_init=.TRUE. 200 ! 201 IF( ln_bt_fw .OR. neuler == 0 ) THEN 202 ll_fw_start =.TRUE. 203 noffset = 0 204 204 ELSE 205 ll_fw_start=.FALSE.205 ll_fw_start =.FALSE. 206 206 ENDIF 207 207 ! 208 208 ! Set averaging weights and cycle length: 209 CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2) 210 ! 209 CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) 211 210 ! 212 211 ENDIF … … 219 218 ! and update depths at T-F points (ht and zhf resp.) at each barotropic time step 220 219 ! 221 IF ( kt == nit000 .OR. lk_vvl) THEN222 IF ( ln_dynvor_een ) THEN!== EEN scheme ==!220 IF( kt == nit000 .OR. .NOT.ln_linssh ) THEN 221 IF( ln_dynvor_een ) THEN !== EEN scheme ==! 223 222 SELECT CASE( nn_een_e3f ) !* ff/e3 at F-point 224 223 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 225 224 DO jj = 1, jpjm1 226 225 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._wp226 zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + & 227 & ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) * 0.25_wp 229 228 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff(ji,jj) / zwz(ji,jj) 230 229 END DO … … 233 232 DO jj = 1, jpjm1 234 233 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 ) ) &234 zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + & 235 & ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) & 237 236 & / ( MAX( 1._wp, tmask(ji ,jj+1, 1) + tmask(ji+1,jj+1, 1) + & 238 237 & tmask(ji ,jj , 1) + tmask(ji+1,jj , 1) ) ) … … 276 275 DO jk = 1, jpkm1 277 276 DO jj = 1, jpjm1 278 zhf(:,jj) = zhf(:,jj) + fse3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)277 zhf(:,jj) = zhf(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 279 278 END DO 280 279 END DO … … 308 307 ! 309 308 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)309 zu_frc(:,:) = zu_frc(:,:) + e3u_n(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 310 zv_frc(:,:) = zv_frc(:,:) + e3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 312 311 END DO 313 312 ! 314 zu_frc(:,:) = zu_frc(:,:) * hur(:,:)315 zv_frc(:,:) = zv_frc(:,:) * hvr(:,:)313 zu_frc(:,:) = zu_frc(:,:) * r1_hu_n(:,:) 314 zv_frc(:,:) = zv_frc(:,:) * r1_hv_n(:,:) 316 315 ! 317 316 ! … … 327 326 ! !* barotropic Coriolis trends (vorticity scheme dependent) 328 327 ! ! -------------------------------------------------------- 329 zwx(:,:) = un_b(:,:) * hu (:,:) * e2u(:,:) ! now fluxes330 zwy(:,:) = vn_b(:,:) * hv (:,:) * e1v(:,:)328 zwx(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:) ! now fluxes 329 zwy(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:) 331 330 ! 332 331 IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN ! energy conserving or mixed scheme … … 373 372 ! !* Right-Hand-Side of the barotropic momentum equation 374 373 ! ! ---------------------------------------------------- 375 IF( lk_vvl ) THEN ! Variable volume : remove surface pressure gradient 376 DO jj = 2, jpjm1 377 DO ji = fs_2, fs_jpim1 ! vector opt. 378 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) * r1_e1u(ji,jj) 379 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) * r1_e2v(ji,jj) 380 END DO 381 END DO 374 IF( .NOT.ln_linssh ) THEN ! Variable volume : remove surface pressure gradient 375 IF( ln_wd ) THEN ! Calculating and applying W/D gravity filters 376 wduflt1(:,:) = 1.0_wp 377 wdvflt1(:,:) = 1.0_wp 378 DO jj = 2, jpjm1 379 DO ji = 2, jpim1 380 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 381 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) & 382 & > rn_wdmin1 + rn_wdmin2 383 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 384 & + rn_wdmin1 + rn_wdmin2 385 IF(ll_tmp1) THEN 386 zcpx(ji,jj) = 1.0_wp 387 ELSEIF(ll_tmp2) THEN 388 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen here 389 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) & 390 & /(sshn(ji+1,jj) - sshn(ji,jj))) 391 ELSE 392 zcpx(ji,jj) = 0._wp 393 wduflt1(ji,jj) = 0.0_wp 394 END IF 395 396 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 397 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) & 398 & > rn_wdmin1 + rn_wdmin2 399 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 400 & + rn_wdmin1 + rn_wdmin2 401 IF(ll_tmp1) THEN 402 zcpy(ji,jj) = 1.0_wp 403 ELSEIF(ll_tmp2) THEN 404 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen here 405 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) & 406 & /(sshn(ji,jj+1) - sshn(ji,jj))) 407 ELSE 408 zcpy(ji,jj) = 0._wp 409 wdvflt1(ji,jj) = 0.0_wp 410 ENDIF 411 412 END DO 413 END DO 414 415 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 416 417 DO jj = 2, jpjm1 418 DO ji = 2, jpim1 419 zu_trd(ji,jj) = ( zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) & 420 & * r1_e1u(ji,jj) ) * zcpx(ji,jj) * wduflt1(ji,jj) 421 zv_trd(ji,jj) = ( zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) & 422 & * r1_e2v(ji,jj) ) * zcpy(ji,jj) * wdvflt1(ji,jj) 423 END DO 424 END DO 425 426 ELSE 427 428 DO jj = 2, jpjm1 429 DO ji = fs_2, fs_jpim1 ! vector opt. 430 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) * r1_e1u(ji,jj) 431 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) * r1_e2v(ji,jj) 432 END DO 433 END DO 434 ENDIF 435 382 436 ENDIF 383 437 384 438 DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend 385 439 DO ji = fs_2, fs_jpim1 386 zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * umask(ji,jj,1)387 zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1)440 zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj) 441 zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj) 388 442 END DO 389 443 END DO … … 411 465 ! 412 466 ! 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(:,:) 467 IF( ln_wd ) THEN 468 zu_frc(:,:) = zu_frc(:,:) + MAX(r1_hu_n(:,:) * bfrua(:,:),-1._wp / rdtbt) * zwx(:,:) 469 zv_frc(:,:) = zv_frc(:,:) + MAX(r1_hv_n(:,:) * bfrva(:,:),-1._wp / rdtbt) * zwy(:,:) 470 ELSE 471 zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * bfrua(:,:) * zwx(:,:) 472 zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * bfrva(:,:) * zwy(:,:) 473 END IF 474 ! 475 ! ! Add top stress contribution from baroclinic velocities: 476 IF (ln_bt_fw) THEN 477 DO jj = 2, jpjm1 478 DO ji = fs_2, fs_jpim1 ! vector opt. 479 iktu = miku(ji,jj) 480 iktv = mikv(ji,jj) 481 zwx(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) ! NOW top baroclinic velocities 482 zwy(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj) 483 END DO 484 END DO 485 ELSE 486 DO jj = 2, jpjm1 487 DO ji = fs_2, fs_jpim1 ! vector opt. 488 iktu = miku(ji,jj) 489 iktv = mikv(ji,jj) 490 zwx(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) ! BEFORE top baroclinic velocities 491 zwy(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj) 492 END DO 493 END DO 494 ENDIF 495 ! 496 ! Note that the "unclipped" top friction parameter is used even with explicit drag 497 zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * tfrua(:,:) * zwx(:,:) 498 zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * tfrva(:,:) * zwy(:,:) 415 499 ! 416 500 IF (ln_bt_fw) THEN ! Add wind forcing 417 zu_frc(:,:) = zu_frc(:,:) + zraur * utau(:,:) * hur(:,:)418 zv_frc(:,:) = zv_frc(:,:) + zraur * vtau(:,:) * hvr(:,:)501 zu_frc(:,:) = zu_frc(:,:) + zraur * utau(:,:) * r1_hu_n(:,:) 502 zv_frc(:,:) = zv_frc(:,:) + zraur * vtau(:,:) * r1_hv_n(:,:) 419 503 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(:,:)504 zu_frc(:,:) = zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * r1_hu_n(:,:) 505 zv_frc(:,:) = zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * r1_hv_n(:,:) 422 506 ENDIF 423 507 ! … … 482 566 vb_e (:,:) = 0._wp 483 567 ENDIF 568 569 IF( ln_wd ) THEN !preserve the positivity of water depth 570 !ssh[b,n,a] should have already been processed for this 571 sshbb_e(:,:) = MAX(sshbb_e(:,:), rn_wdmin1 - bathy(:,:)) 572 sshb_e(:,:) = MAX(sshb_e(:,:) , rn_wdmin1 - bathy(:,:)) 573 ENDIF 484 574 ! 485 575 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(:,:)576 sshn_e(:,:) = sshn(:,:) 577 un_e (:,:) = un_b(:,:) 578 vn_e (:,:) = vn_b(:,:) 579 ! 580 hu_e (:,:) = hu_n(:,:) 581 hv_e (:,:) = hv_n(:,:) 582 hur_e (:,:) = r1_hu_n(:,:) 583 hvr_e (:,:) = r1_hv_n(:,:) 494 584 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(:,:)585 sshn_e(:,:) = sshb(:,:) 586 un_e (:,:) = ub_b(:,:) 587 vn_e (:,:) = vb_b(:,:) 588 ! 589 hu_e (:,:) = hu_b(:,:) 590 hv_e (:,:) = hv_b(:,:) 591 hur_e (:,:) = r1_hu_b(:,:) 592 hvr_e (:,:) = r1_hv_b(:,:) 503 593 ENDIF 504 594 ! … … 518 608 ! Update only tidal forcing at open boundaries 519 609 #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)610 IF( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 ) 611 IF( ln_tide_pot .AND. lk_tide ) CALL upd_tide ( kt, kit=jn, time_offset= noffset ) 522 612 #endif 523 613 ! … … 537 627 va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:) 538 628 539 IF( lk_vvl ) THEN!* Update ocean depth (variable volume case only)629 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 540 630 ! ! ------------------ 541 631 ! Extrapolate Sea Level at step jit+0.5: … … 544 634 DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points 545 635 DO ji = 2, fs_jpim1 ! Vector opt. 546 zwx(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) &636 zwx(ji,jj) = z1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 547 637 & * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 548 638 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) 549 zwy(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) &639 zwy(ji,jj) = z1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 550 640 & * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 551 641 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) … … 556 646 zhup2_e (:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points 557 647 zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:) 648 IF( ln_wd ) THEN 649 zhup2_e(:,:) = MAX(zhup2_e (:,:), rn_wdmin1) 650 zhvp2_e(:,:) = MAX(zhvp2_e (:,:), rn_wdmin1) 651 END IF 558 652 ELSE 559 zhup2_e (:,:) = hu (:,:)560 zhvp2_e (:,:) = hv (:,:)653 zhup2_e (:,:) = hu_n(:,:) 654 zhvp2_e (:,:) = hv_n(:,:) 561 655 ENDIF 562 656 ! !* after ssh … … 569 663 ! 570 664 #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 665 ! Set fluxes during predictor step to ensure volume conservation 666 IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN 574 667 IF((nbondi == -1).OR.(nbondi == 2)) THEN 575 668 DO jj=1,jpj … … 594 687 ENDIF 595 688 #endif 689 IF( ln_wd ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 596 690 ! 597 691 ! Sum over sub-time-steps to compute advective velocities … … 607 701 END DO 608 702 END DO 609 ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1) 703 ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * ssmask(:,:) 704 IF( ln_wd ) ssha_e(:,:) = MAX(ssha_e(:,:), rn_wdmin1 - bathy(:,:)) 610 705 CALL lbc_lnk( ssha_e, 'T', 1._wp ) 611 706 612 707 #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 )708 ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T) 709 IF( lk_bdy ) CALL bdy_ssh( ssha_e ) 615 710 #endif 616 711 #if defined key_agrif 617 IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )712 IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn ) 618 713 #endif 619 714 ! 620 715 ! Sea Surface Height at u-,v-points (vvl case only) 621 IF ( lk_vvl) THEN716 IF( .NOT.ln_linssh ) THEN 622 717 DO jj = 2, jpjm1 623 718 DO ji = 2, jpim1 ! NO Vector Opt. 624 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj)&625 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &626 & + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) )627 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj)&628 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) &629 & + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) )719 zsshu_a(ji,jj) = z1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 720 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 721 & + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) ) 722 zsshv_a(ji,jj) = z1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 723 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 724 & + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) ) 630 725 END DO 631 726 END DO … … 651 746 za3=0.013_wp ! za3 = eps 652 747 ENDIF 653 748 ! 654 749 zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) & 655 750 & + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 656 751 IF( ln_wd ) THEN ! Calculating and applying W/D gravity filters 752 wduflt1(:,:) = 1._wp 753 wdvflt1(:,:) = 1._wp 754 DO jj = 2, jpjm1 755 DO ji = 2, jpim1 756 ll_tmp1 = MIN( zsshp2_e(ji,jj), zsshp2_e(ji+1,jj) ) > MAX( -bathy(ji,jj), -bathy(ji+1,jj) ) & 757 & .AND. MAX( zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji+1,jj) + bathy(ji+1,jj) ) & 758 & > rn_wdmin1 + rn_wdmin2 759 ll_tmp2 = MAX( zsshp2_e(ji,jj), zsshp2_e(ji+1,jj) ) > MAX( -bathy(ji,jj), -bathy(ji+1,jj) ) & 760 & + rn_wdmin1 + rn_wdmin2 761 IF(ll_tmp1) THEN 762 zcpx(ji,jj) = 1._wp 763 ELSE IF(ll_tmp2) THEN 764 ! no worries about zsshp2_e(ji+1,jj)-zsshp2_e(ji,jj) = 0, it won't happen here 765 zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + bathy(ji+1,jj) - zsshp2_e(ji,jj) - bathy(ji,jj)) & 766 & / (zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj)) ) 767 ELSE 768 zcpx(ji,jj) = 0._wp 769 wduflt1(ji,jj) = 0._wp 770 END IF 771 772 ll_tmp1 = MIN( zsshp2_e(ji,jj), zsshp2_e(ji,jj+1) ) > MAX( -bathy(ji,jj), -bathy(ji,jj+1) ) & 773 & .AND. MAX( zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji,jj+1) + bathy(ji,jj+1) ) & 774 & > rn_wdmin1 + rn_wdmin2 775 ll_tmp2 = MAX( zsshp2_e(ji,jj), zsshp2_e(ji,jj+1) ) > MAX( -bathy(ji,jj), -bathy(ji,jj+1) ) & 776 & + rn_wdmin1 + rn_wdmin2 777 IF(ll_tmp1) THEN 778 zcpy(ji,jj) = 1._wp 779 ELSE IF(ll_tmp2) THEN 780 ! no worries about zsshp2_e(ji,jj+1)-zsshp2_e(ji,jj) = 0, it won't happen here 781 zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + bathy(ji,jj+1) - zsshp2_e(ji,jj) - bathy(ji,jj)) & 782 & / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj)) ) 783 ELSE 784 zcpy(ji,jj) = 0._wp 785 wdvflt1(ji,jj) = 0._wp 786 END IF 787 END DO 788 END DO 789 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 790 ENDIF 657 791 ! 658 792 ! Compute associated depths at U and V points: 659 IF ( lk_vvl.AND.(.NOT.ln_dynadv_vec) ) THEN793 IF( .NOT.ln_linssh .AND. .NOT.ln_dynadv_vec ) THEN !* Vector form 660 794 ! 661 795 DO jj = 2, jpjm1 662 796 DO ji = 2, jpim1 663 zx1 = z1_2 * umask(ji ,jj,1) * r1_e1e2u(ji ,jj) &797 zx1 = z1_2 * ssumask(ji ,jj) * r1_e1e2u(ji ,jj) & 664 798 & * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj) & 665 799 & + e1e2t(ji+1,jj ) * zsshp2_e(ji+1,jj ) ) 666 zy1 = z1_2 * vmask(ji ,jj,1) * r1_e1e2v(ji ,jj ) &800 zy1 = z1_2 * ssvmask(ji ,jj) * r1_e1e2v(ji ,jj ) & 667 801 & * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj ) & 668 802 & + e1e2t(ji ,jj+1) * zsshp2_e(ji ,jj+1) ) … … 671 805 END DO 672 806 END DO 807 808 IF( ln_wd ) THEN 809 zhust_e(:,:) = MAX(zhust_e (:,:), rn_wdmin1 ) 810 zhvst_e(:,:) = MAX(zhvst_e (:,:), rn_wdmin1 ) 811 END IF 812 673 813 ENDIF 674 814 ! 675 815 ! Add Coriolis trend: 676 ! zwz array below or triads normally depend on sea level with key_vvland should be updated816 ! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated 677 817 ! at each time step. We however keep them constant here for optimization. 678 818 ! Recall that zwx and zwy arrays hold fluxes at this stage: … … 680 820 ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 681 821 ! 682 IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN 822 IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN !== energy conserving or mixed scheme ==! 683 823 DO jj = 2, jpjm1 684 824 DO ji = fs_2, fs_jpim1 ! vector opt. … … 692 832 END DO 693 833 ! 694 ELSEIF ( ln_dynvor_ens ) THEN 834 ELSEIF ( ln_dynvor_ens ) THEN !== enstrophy conserving scheme ==! 695 835 DO jj = 2, jpjm1 696 836 DO ji = fs_2, fs_jpim1 ! vector opt. … … 704 844 END DO 705 845 ! 706 ELSEIF ( ln_dynvor_een ) THEN !== energy and enstrophy conserving scheme ==!846 ELSEIF ( ln_dynvor_een ) THEN !== energy and enstrophy conserving scheme ==! 707 847 DO jj = 2, jpjm1 708 848 DO ji = fs_2, fs_jpim1 ! vector opt. … … 736 876 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 737 877 ! 878 ! Add top stresses: 879 zu_trd(:,:) = zu_trd(:,:) + tfrua(:,:) * un_e(:,:) * hur_e(:,:) 880 zv_trd(:,:) = zv_trd(:,:) + tfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 881 ! 738 882 ! Surface pressure trend: 739 DO jj = 2, jpjm1 740 DO ji = fs_2, fs_jpim1 ! vector opt. 741 ! Add surface pressure gradient 742 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 743 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 744 zwx(ji,jj) = zu_spg 745 zwy(ji,jj) = zv_spg 746 END DO 747 END DO 883 884 IF( ln_wd ) THEN 885 DO jj = 2, jpjm1 886 DO ji = 2, jpim1 887 ! Add surface pressure gradient 888 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 889 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 890 zwx(ji,jj) = zu_spg * zcpx(ji,jj) 891 zwy(ji,jj) = zv_spg * zcpy(ji,jj) 892 END DO 893 END DO 894 ELSE 895 DO jj = 2, jpjm1 896 DO ji = fs_2, fs_jpim1 ! vector opt. 897 ! Add surface pressure gradient 898 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 899 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 900 zwx(ji,jj) = zu_spg 901 zwy(ji,jj) = zv_spg 902 END DO 903 END DO 904 END IF 905 748 906 ! 749 907 ! Set next velocities: 750 IF( ln_dynadv_vec .OR. (.NOT. lk_vvl) ) THEN !Vector form908 IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form 751 909 DO jj = 2, jpjm1 752 910 DO ji = fs_2, fs_jpim1 ! vector opt. … … 755 913 & + zu_trd(ji,jj) & 756 914 & + zu_frc(ji,jj) ) & 757 & ) * umask(ji,jj,1)915 & ) * ssumask(ji,jj) 758 916 759 917 va_e(ji,jj) = ( vn_e(ji,jj) & … … 761 919 & + zv_trd(ji,jj) & 762 920 & + zv_frc(ji,jj) ) & 763 & ) * vmask(ji,jj,1)764 END DO 765 END DO 766 767 ELSE !Flux form921 & ) * ssvmask(ji,jj) 922 END DO 923 END DO 924 ! 925 ELSE !* Flux form 768 926 DO jj = 2, jpjm1 769 927 DO ji = fs_2, fs_jpim1 ! vector opt. 770 928 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)) 929 IF( ln_wd ) THEN 930 zhura = MAX(hu_0(ji,jj) + zsshu_a(ji,jj), rn_wdmin1) 931 zhvra = MAX(hv_0(ji,jj) + zsshv_a(ji,jj), rn_wdmin1) 932 ELSE 933 zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 934 zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 935 END IF 936 zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj)) 937 zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj)) 773 938 774 939 ua_e(ji,jj) = ( hu_e(ji,jj) * un_e(ji,jj) & 775 940 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 776 941 & + zhup2_e(ji,jj) * zu_trd(ji,jj) & 777 & + hu(ji,jj) * zu_frc(ji,jj) ) &942 & + hu_n(ji,jj) * zu_frc(ji,jj) ) & 778 943 & ) * zhura 779 944 … … 781 946 & + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) & 782 947 & + zhvp2_e(ji,jj) * zv_trd(ji,jj) & 783 & + hv(ji,jj) * zv_frc(ji,jj) ) &948 & + hv_n(ji,jj) * zv_frc(ji,jj) ) & 784 949 & ) * zhvra 785 950 END DO … … 787 952 ENDIF 788 953 ! 789 IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only) 790 ! ! ---------------------------------------------- 791 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 792 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 793 hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) ) 794 hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) ) 954 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 955 IF( ln_wd ) THEN 956 hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1) 957 hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1) 958 ELSE 959 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 960 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 961 END IF 962 hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) ) 963 hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) ) 795 964 ! 796 965 ENDIF 797 ! !* domain lateral boundary 798 ! ! ----------------------- 799 ! 966 ! !* domain lateral boundary 800 967 CALL lbc_lnk_multi( ua_e, 'U', -1._wp, va_e , 'V', -1._wp ) 801 968 ! 802 969 #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 )970 ! ! open boundaries 971 IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e ) 805 972 #endif 806 973 #if defined key_agrif … … 824 991 ! ! ---------------------- 825 992 za1 = wgtbtp1(jn) 826 IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN ! Sum velocities993 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! Sum velocities 827 994 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) 828 995 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) 829 ELSE 996 ELSE ! Sum transports 830 997 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) 831 998 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) … … 843 1010 zwx(:,:) = un_adv(:,:) 844 1011 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(:,:)1012 IF( ( kt == nit000 .AND. neuler==0 ) .OR. .NOT.ln_bt_fw ) THEN 1013 un_adv(:,:) = zwx(:,:) * r1_hu_n(:,:) 1014 vn_adv(:,:) = zwy(:,:) * r1_hv_n(:,:) 848 1015 ELSE 849 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:) ) * hur(:,:)850 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zwy(:,:) ) * hvr(:,:)1016 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:) ) * r1_hu_n(:,:) 1017 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zwy(:,:) ) * r1_hv_n(:,:) 851 1018 END IF 852 1019 853 IF (ln_bt_fw) THEN ! Save integrated transport for next computation1020 IF( ln_bt_fw ) THEN ! Save integrated transport for next computation 854 1021 ub2_b(:,:) = zwx(:,:) 855 1022 vb2_b(:,:) = zwy(:,:) … … 857 1024 ! 858 1025 ! Update barotropic trend: 859 IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN1026 IF( ln_dynadv_vec .OR. ln_linssh ) THEN 860 1027 DO jk=1,jpkm1 861 1028 ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b … … 877 1044 ! 878 1045 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_b1046 ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b 1047 va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b 881 1048 END DO 882 1049 ! Save barotropic velocities not transport: 883 ua_b (:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - umask(:,:,1) )884 va_b (:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - vmask(:,:,1) )1050 ua_b(:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) ) 1051 va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 885 1052 ENDIF 886 1053 ! 887 1054 DO jk = 1, jpkm1 888 1055 ! Correct velocities: 889 un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) ) *umask(:,:,jk)890 vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) ) *vmask(:,:,jk)1056 un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) ) * umask(:,:,jk) 1057 vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) ) * vmask(:,:,jk) 891 1058 ! 892 1059 END DO 1060 ! 1061 CALL iom_put( "ubar", un_adv(:,:) ) ! barotropic i-current 1062 CALL iom_put( "vbar", vn_adv(:,:) ) ! barotropic i-current 893 1063 ! 894 1064 #if defined key_agrif … … 896 1066 ! (used to update coarse grid transports at next time step) 897 1067 ! 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. e01068 IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN 1069 IF( Agrif_NbStepint() == 0 ) THEN 1070 ub2_i_b(:,:) = 0._wp 1071 vb2_i_b(:,:) = 0._wp 902 1072 END IF 903 1073 ! … … 906 1076 vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:) 907 1077 ENDIF 908 !909 !910 1078 #endif 911 !912 1079 ! !* 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 ) 921 ! 1080 IF( lrst_oce .AND.ln_bt_fw ) CALL ts_rst( kt, 'WRITE' ) 1081 ! 1082 CALL wrk_dealloc( jpi,jpj, zsshp2_e, zhdiv ) 1083 CALL wrk_dealloc( jpi,jpj, zu_trd, zv_trd ) 1084 CALL wrk_dealloc( jpi,jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc ) 1085 CALL wrk_dealloc( jpi,jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e ) 1086 CALL wrk_dealloc( jpi,jpj, zsshu_a, zsshv_a ) 1087 CALL wrk_dealloc( jpi,jpj, zhf ) 1088 IF( ln_wd ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 ) 1089 ! 1090 IF ( ln_diatmb ) THEN 1091 CALL iom_put( "baro_u" , un_b*umask(:,:,1)+zmdi*(1-umask(:,:,1 ) ) ) ! Barotropic U Velocity 1092 CALL iom_put( "baro_v" , vn_b*vmask(:,:,1)+zmdi*(1-vmask(:,:,1 ) ) ) ! Barotropic V Velocity 1093 ENDIF 922 1094 IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_ts') 923 1095 ! 924 1096 END SUBROUTINE dyn_spg_ts 1097 925 1098 926 1099 SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2) … … 1001 1174 END SUBROUTINE ts_wgt 1002 1175 1176 1003 1177 SUBROUTINE ts_rst( kt, cdrw ) 1004 1178 !!--------------------------------------------------------------------- … … 1054 1228 END SUBROUTINE ts_rst 1055 1229 1056 SUBROUTINE dyn_spg_ts_init( kt ) 1230 1231 SUBROUTINE dyn_spg_ts_init 1057 1232 !!--------------------------------------------------------------------- 1058 1233 !! *** ROUTINE dyn_spg_ts_init *** … … 1060 1235 !! ** Purpose : Set time splitting options 1061 1236 !!---------------------------------------------------------------------- 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 !! 1237 INTEGER :: ji ,jj ! dummy loop indices 1238 REAL(wp) :: zxr2, zyr2, zcmax ! local scalar 1239 REAL(wp), POINTER, DIMENSION(:,:) :: zcu 1068 1240 !!---------------------------------------------------------------------- 1069 1241 ! 1070 1242 ! Max courant number for ext. grav. waves 1071 1243 ! 1072 CALL wrk_alloc( jpi, jpj,zcu )1244 CALL wrk_alloc( jpi,jpj, zcu ) 1073 1245 ! 1074 1246 DO jj = 1, jpj … … 1084 1256 1085 1257 ! 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)1258 IF( ln_bt_auto ) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax) 1087 1259 1088 1260 rdtbt = rdt / REAL( nn_baro , wp ) … … 1114 1286 #if defined key_agrif 1115 1287 ! 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.' )1288 IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' ) 1117 1289 #endif 1118 1290 ! 1119 1291 IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt 1120 1292 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 1293 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac' 1294 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_baro' 1295 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = 2*nn_baro' 1296 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' ) 1125 1297 END SELECT 1126 1298 ! … … 1130 1302 IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax 1131 1303 ! 1132 IF ((.NOT.ln_bt_av).AND.(.NOT.ln_bt_fw)) THEN1304 IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN 1133 1305 CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' ) 1134 1306 ENDIF 1135 IF 1307 IF( zcmax>0.9_wp ) THEN 1136 1308 CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' ) 1137 1309 ENDIF 1138 1310 ! 1139 CALL wrk_dealloc( jpi, jpj,zcu )1311 CALL wrk_dealloc( jpi,jpj, zcu ) 1140 1312 ! 1141 1313 END SUBROUTINE dyn_spg_ts_init
Note: See TracChangeset
for help on using the changeset viewer.
