Changeset 4241
- Timestamp:
- 2013-11-19T09:04:21+01:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r4194 r4241 9 9 !! 3.3 ! 2010-09 (D. Storkey, E. O'Dea) update for BDY for Shelf configurations 10 10 !! 3.3 ! 2011-03 (R. Benshila, R. Hordoir, P. Oddo) update calculation of ub_b 11 !! 3.5 ! 2013-07 (J. Chanut) Switch to Forward-backward time stepping 12 !! 3.6 ! 2013-11 (A. Coward) Update for z-tilde compatibility 11 13 !!--------------------------------------------------------------------- 12 14 #if defined key_dynspg_ts || defined key_esopa … … 16 18 !! dyn_spg_ts : compute surface pressure gradient trend using a time- 17 19 !! splitting scheme and add to the general trend 18 !! ts_rst : read/write the time-splitting restart fields in the ocean restart file19 20 !!---------------------------------------------------------------------- 20 21 USE oce ! ocean dynamics and tracers … … 24 25 USE phycst ! physical constants 25 26 USE domvvl ! variable volume 26 USE zdfbfr ! bottom friction27 27 USE dynvor ! vorticity term 28 USE dynadv ! advection29 USE obc_oce ! Lateral open boundary condition30 USE obc_par ! open boundary condition parameters31 USE obcdta ! open boundary condition data32 USE obcfla ! Flather open boundary condition33 28 USE bdy_par ! for lk_bdy 34 29 USE bdy_oce ! Lateral open boundary condition 35 USE bdy dta! open boundary condition data30 USE bdytides ! open boundary condition data 36 31 USE bdydyn2d ! open boundary conditions on barotropic variables 37 USE bdytides ! 38 USE sbctide 39 USE updtide 32 USE sbctide ! tides 33 USE updtide ! tide potential 40 34 USE lib_mpp ! distributed memory computing library 41 35 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 43 37 USE in_out_manager ! I/O manager 44 38 USE iom ! IOM library 39 USE restart ! only for lrst_oce 45 40 USE zdf_oce ! Vertical diffusion 46 41 USE wrk_nemo ! Memory Allocation 47 USE timing ! Timing 42 USE timing ! Timing 43 USE sbcapr ! surface boundary condition: atmospheric pressure 44 USE dynadv, ONLY: ln_dynadv_vec 48 45 49 46 50 47 IMPLICIT NONE 51 48 PRIVATE 49 50 PUBLIC dyn_spg_ts ! routine called in dynspg.F90 51 PUBLIC dyn_spg_ts_alloc ! " " " " 52 PUBLIC dyn_spg_ts_init ! " " " " 52 53 53 54 ! Potential namelist parameters below to be read in dyn_spg_ts_init … … 65 66 wgtbtp1, & ! Primary weights used for time filtering of barotropic variables 66 67 wgtbtp2 ! Secondary weights used for time filtering of barotropic variables 67 !68 PUBLIC dyn_spg_ts ! routine called by step.F9069 PUBLIC ts_rst ! routine called by istate.F9070 PUBLIC dyn_spg_ts_alloc ! routine called by dynspg.F9071 PUBLIC dyn_spg_ts_init ! " " " "72 73 68 74 69 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz ! ff/h at F points … … 96 91 # include "vectopt_loop_substitute.h90" 97 92 !!---------------------------------------------------------------------- 98 !! NEMO/OPA 4.0 , NEMO Consortium (2011)99 !! $Id $93 !! NEMO/OPA 3.5 , NEMO Consortium (2013) 94 !! $Id: dynspg_ts.F90 100 95 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 101 96 !!---------------------------------------------------------------------- … … 132 127 END FUNCTION dyn_spg_ts_alloc 133 128 134 135 129 SUBROUTINE dyn_spg_ts( kt ) 136 130 !!---------------------------------------------------------------------- 137 !! *** routine dyn_spg_ts ***138 131 !! 139 !! ** Purpose : Compute the now trend due to the surface pressure 140 !! gradient in case of free surface formulation with time-splitting. 141 !! Add it to the general trend of momentum equation. 132 !! ** Purpose : 133 !! -Compute the now trend due to the explicit time stepping 134 !! of the quasi-linear barotropic system. Barotropic variables are 135 !! advanced from internal time steps "n" to "n+1" (if ln_bt_cen=F) 136 !! or from "n-1" to "n+1" time steps (if ln_bt_cen=T) with a 137 !! generalized forward-backward (see ref. below) time stepping. 138 !! -Update the free surface at step "n+1" (ssha, zsshu_a, zsshv_a). 139 !! -Compute barotropic advective velocities at step "n" to be used 140 !! to advect tracers latter on. These are compliant with discrete 141 !! continuity equation taken at the baroclinic time steps, thus 142 !! ensuring tracers conservation. 142 143 !! 143 !! ** Method : Free surface formulation with time-splitting 144 !! -1- Save the vertically integrated trend. This general trend is 145 !! held constant over the barotropic integration. 146 !! The Coriolis force is removed from the general trend as the 147 !! surface gradient and the Coriolis force are updated within 148 !! the barotropic integration. 149 !! -2- Barotropic loop : updates of sea surface height (ssha_e) and 150 !! barotropic velocity (ua_e and va_e) through barotropic 151 !! momentum and continuity integration. Barotropic former 152 !! variables are time averaging over the full barotropic cycle 153 !! (= 2 * baroclinic time step) and saved in uX_b 154 !! and vX_b (X specifying after, now or before). 155 !! -3- The new general trend becomes : 156 !! ua = ua - sum_k(ua)/H + ( un_b - ub_b ) 144 !! ** Method : 157 145 !! 158 !! ** Action : - Update (ua,va) with the surf. pressure gradient trend 146 !! ** Action : - Update barotropic velocities: ua_b, va_b 147 !! - Update trend (ua,va) with barotropic component 148 !! - Update ssha, zsshu_a, zsshv_a 149 !! - Update barotropic advective velocity at kt=now 159 150 !! 160 !! References : Griffies et al., (2003): A technical guide to MOM4. NOAA/GFDL 151 !! References : Shchepetkin, A.F. and J.C. McWilliams, 2005: 152 !! The regional oceanic modeling system (ROMS): 153 !! a split-explicit, free-surface, 154 !! topography-following-coordinate oceanic model. 155 !! Ocean Modelling, 9, 347-404. 161 156 !!--------------------------------------------------------------------- 162 157 ! 163 158 INTEGER, INTENT(in) :: kt ! ocean time-step index 164 159 ! 165 LOGICAL :: ll_fw_start ! if true, forward integration 166 LOGICAL :: ll_init ! if true, special startup of 2d equations 167 INTEGER :: ji, jj, jk, jn ! dummy loop indices 168 ! INTEGER :: icycle ! local scalar 169 INTEGER :: ioffset ! local scalar 170 INTEGER :: ikbu, ikbv ! local scalar 171 REAL(wp) :: zraur, zcoef, z2dt_e, z1_2dt_b, z2dt_bf ! local scalars 172 REAL(wp) :: zx1, zy1, zx2, zy2 ! - - 173 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 174 REAL(wp) :: za0, za1, za2, za3 ! - - 175 REAL(wp) :: zu_spg, zu_cor, zu_sld, zu_asp ! - - 176 REAL(wp) :: zv_spg, zv_cor, zv_sld, zv_asp ! - - 177 REAL(wp) :: ua_btm, va_btm ! - - 178 ! 179 REAL(wp), POINTER, DIMENSION(:,:) :: zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv 180 REAL(wp), POINTER, DIMENSION(:,:) :: zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e 181 REAL(wp), POINTER, DIMENSION(:,:) :: zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum 182 REAL(wp), POINTER, DIMENSION(:,:) :: zhu_b, zhv_b 183 REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zsshp2_e 184 REAL(wp), POINTER, DIMENSION(:,:) :: zhust_e, zhvst_e 160 LOGICAL :: ll_fw_start ! if true, forward integration 161 LOGICAL :: ll_init ! if true, special startup of 2d equations 162 INTEGER :: ji, jj, jk, jn ! dummy loop indices 163 INTEGER :: ikbu, ikbv, noffset ! local integers 164 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars 165 REAL(wp) :: zx1, zy1, zx2, zy2 ! - - 166 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 ! - - 167 REAL(wp) :: zu_spg, zv_spg ! - - 168 REAL(wp) :: zhura, zhvra ! - - 169 REAL(wp) :: za0, za1, za2, za3 ! - - 170 ! 171 REAL(wp), POINTER, DIMENSION(:,:) :: zun_e, zvn_e, zsshp2_e 172 REAL(wp), POINTER, DIMENSION(:,:) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc 173 REAL(wp), POINTER, DIMENSION(:,:) :: zu_sum, zv_sum, zwx, zwy, zhdiv 174 REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e 175 REAL(wp), POINTER, DIMENSION(:,:) :: zhur_b, zhvr_b 176 REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a 177 REAL(wp), POINTER, DIMENSION(:,:) :: zht, zhf 185 178 !!---------------------------------------------------------------------- 186 179 ! 187 180 IF( nn_timing == 1 ) CALL timing_start('dyn_spg_ts') 188 181 ! 189 CALL wrk_alloc( jpi, jpj, zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv ) 190 CALL wrk_alloc( jpi, jpj, zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e ) 191 CALL wrk_alloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum ) 192 CALL wrk_alloc( jpi, jpj, zhu_b, zhv_b ) 193 CALL wrk_alloc( jpi, jpj, zhup2_e, zhvp2_e, zsshp2_e ) 194 CALL wrk_alloc( jpi, jpj, zhust_e, zhvst_e ) 195 ! 196 197 ! !* Local constant initialization 198 z1_12 = 1._wp / 12._wp 199 z1_8 = 0.125_wp ! coefficients for vorticity estimates 200 z1_4 = 0.25_wp 201 z1_2 = 0.5_wp 202 zraur = 1._wp / rau0 ! 1 / volumetric mass 203 ! 204 zhdiv(:,:) = 0._wp ! barotropic divergence 205 zu_sld = 0._wp ; zu_asp = 0._wp ! tides trends (lk_tide=F) 206 zv_sld = 0._wp ; zv_asp = 0._wp 207 208 IF( kt == nit000 .AND. neuler == 0) THEN ! for implicit bottom friction 209 z2dt_bf = rdt ! baroclinic time step (euler timestep) 210 ELSE 211 z2dt_bf = 2.0_wp * rdt ! baroclinic time step 212 ENDIF 213 z1_2dt_b = 1.0_wp / z2dt_bf ! reciprocal of baroclinic time step 214 ! 215 ll_init = ln_bt_av ! if no time averaging, then no specific restart 182 ! !* Allocate temporay arrays 183 CALL wrk_alloc( jpi, jpj, zsshp2_e, zhdiv ) 184 CALL wrk_alloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e ) 185 CALL wrk_alloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc) 186 CALL wrk_alloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e) 187 CALL wrk_alloc( jpi, jpj, zhur_b, zhvr_b ) 188 CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a ) 189 CALL wrk_alloc( jpi, jpj, zht, zhf ) 190 ! 191 ! !* Local constant initialization 192 z1_12 = 1._wp / 12._wp 193 z1_8 = 0.125_wp 194 z1_4 = 0.25_wp 195 z1_2 = 0.5_wp 196 zraur = 1._wp / rau0 197 ! 198 IF( kt == nit000 .AND. neuler == 0 ) THEN ! reciprocal of baroclinic time step 199 z2dt_bf = rdt 200 ELSE 201 z2dt_bf = 2.0_wp * rdt 202 ENDIF 203 z1_2dt_b = 1.0_wp / z2dt_bf 204 ! 205 ll_init = ln_bt_av ! if no time averaging, then no specific restart 216 206 ll_fw_start = .FALSE. 217 207 ! 218 ! time offset in steps for bdy data update219 IF (.NOT.ln_bt_fw) THEN ; ioffset=-2*nn_baro ; ELSE ; ioffset = 0 ; ENDIF220 ! 221 IF( kt == nit000 ) THEN 208 ! time offset in steps for bdy data update 209 IF (.NOT.ln_bt_fw) THEN ; noffset=-2*nn_baro ; ELSE ; noffset = 0 ; ENDIF 210 ! 211 IF( kt == nit000 ) THEN !* initialisation 222 212 ! 223 213 IF(lwp) WRITE(numout,*) … … 230 220 IF (ln_bt_fw.OR.(neuler==0)) THEN 231 221 ll_fw_start=.TRUE. 232 ioffset = 0222 noffset = 0 233 223 ELSE 234 224 ll_fw_start=.FALSE. … … 238 228 CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2) 239 229 ! 240 IF ((neuler/=0).AND.(ln_bt_fw)) CALL ts_rst( nit000, 'READ' ) ! read or initialize the following fields: un_b, vn_b 241 ! 242 ua_e (:,:) = un_b (:,:) 243 va_e (:,:) = vn_b (:,:) 244 hu_e (:,:) = hu (:,:) 245 hv_e (:,:) = hv (:,:) 246 hur_e (:,:) = hur (:,:) 247 hvr_e (:,:) = hvr (:,:) 248 IF( ln_dynvor_een ) THEN 249 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 230 IF ((neuler/=0).AND.(ln_bt_fw)) CALL ts_rst( nit000, 'READ' ) 231 ! 232 ENDIF 233 ! 234 ! Set arrays to remove/compute coriolis trend. 235 ! Do it once at kt=nit000 if volume is fixed, else at each long time step. 236 ! Note that these arrays are also used during barotropic loop. These are however frozen 237 ! although they should be updated in variable volume case. Not a big approximation. 238 ! To remove this approximation, copy lines below inside barotropic loop 239 ! and update depths at T-F points (ht and hf resp.) at each barotropic time step 240 ! 241 IF ( kt == nit000 .OR. lk_vvl ) THEN 242 IF ( ln_dynvor_een ) THEN 243 zht(:,:) = 0. 244 DO jk = 1, jpk 245 zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 246 END DO 247 DO jj = 1, jpjm1 248 DO ji = 1, jpim1 249 zwz(ji,jj) = ( zht(ji ,jj+1) + zht(ji+1,jj+1) + & 250 & zht(ji ,jj ) + zht(ji+1,jj ) ) & 251 & / ( MAX( 1.0_wp, tmask(ji ,jj+1, 1) + tmask(ji+1,jj+1, 1) + & 252 & tmask(ji ,jj , 1) + tmask(ji+1,jj , 1) ) ) 253 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zwz(ji,jj) 254 END DO 255 END DO 256 CALL lbc_lnk( zwz, 'F', 1._wp ) 257 zwz(:,:) = ff(:,:) * zwz(:,:) 258 259 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 250 260 DO jj = 2, jpj 251 261 DO ji = fs_2, jpi ! vector opt. 252 ftne(ji,jj) = ( ff(ji-1,jj ) + ff(ji ,jj ) + ff(ji ,jj-1) ) / 3._wp 253 ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj ) + ff(ji ,jj ) ) / 3._wp 254 ftse(ji,jj) = ( ff(ji ,jj ) + ff(ji ,jj-1) + ff(ji-1,jj-1) ) / 3._wp 255 ftsw(ji,jj) = ( ff(ji ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj ) ) / 3._wp 256 END DO 257 END DO 258 ENDIF 259 ! 262 ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 263 ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) 264 ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) 265 ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) 266 END DO 267 END DO 268 ELSE 269 zwz(:,:) = 0._wp 270 zht(:,:) = 0. 271 IF ( .not. ln_sco ) THEN 272 ! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level 273 ! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth 274 ! ENDIF 275 ! zht(:,:) = gdepw_0(:,:,jk+1) 276 ELSE 277 zht(:,:) = hbatf(:,:) 278 END IF 279 280 DO jj = 1, jpjm1 281 zht(:,jj) = zht(:,jj)*(1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 282 END DO 283 284 DO jk = 1, jpkm1 285 DO jj = 1, jpjm1 286 zht(:,jj) = zht(:,jj) + fse3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 287 END DO 288 END DO 289 CALL lbc_lnk( zht, 'F', 1._wp ) 290 ! JC: TBC. hf should be greater than 0 291 DO jj = 1, jpj 292 DO ji = 1, jpi 293 IF( zht(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zht(ji,jj) ! zht is actually hf here but it saves an array 294 END DO 295 END DO 296 zwz(:,:) = ff(:,:) * zwz(:,:) 297 ENDIF 260 298 ENDIF 261 299 ! … … 267 305 ENDIF 268 306 307 ! before inverse water column height at u- and v- points 308 IF( lk_vvl ) THEN 309 zhur_b(:,:) = 0. 310 zhvr_b(:,:) = 0. 311 DO jk = 1, jpk 312 zhur_b(:,:) = zhur_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk) 313 zhvr_b(:,:) = zhvr_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 314 END DO 315 zhur_b(:,:) = umask(:,:,1) / ( zhur_b(:,:) + 1. - umask(:,:,1) ) 316 zhvr_b(:,:) = vmask(:,:,1) / ( zhvr_b(:,:) + 1. - vmask(:,:,1) ) 317 ELSE 318 zhur_b(:,:) = hur(:,:) 319 zhvr_b(:,:) = hvr(:,:) 320 ENDIF 321 269 322 ! ----------------------------------------------------------------------------- 270 323 ! Phase 1 : Coupling between general trend and barotropic estimates (1st step) 271 324 ! ----------------------------------------------------------------------------- 272 325 ! 326 ! Some vertical sums (at now and before time steps) below could be suppressed 327 ! if one swap barotropic arrays somewhere 328 ! 273 329 ! !* e3*d/dt(Ua), e3*Ub, e3*Vn (Vertically integrated) 274 ! ! -------------------------- 275 zu a(:,:) = 0._wp ; zun(:,:) = 0._wp ; ub_b(:,:) = 0._wp276 zv a(:,:) = 0._wp ; zvn(:,:) = 0._wp ; vb_b(:,:) = 0._wp330 ! ! -------------------------------------------------- 331 zu_frc(:,:) = 0._wp ; ub_b(:,:) = 0._wp ; un_b(:,:) = 0._wp 332 zv_frc(:,:) = 0._wp ; vb_b(:,:) = 0._wp ; vn_b(:,:) = 0._wp 277 333 ! 278 334 DO jk = 1, jpkm1 … … 284 340 DO ji = 1, jpi 285 341 #endif 286 ! ! now trend287 zu a(ji,jj) = zua(ji,jj) + fse3u_n(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)288 zv a(ji,jj) = zva(ji,jj) + fse3v_n(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)289 ! ! now velocity290 zun(ji,jj) = zun(ji,jj) + fse3u_n(ji,jj,jk) * un(ji,jj,jk)291 zvn(ji,jj) = zvn(ji,jj) + fse3v_n(ji,jj,jk) * vn(ji,jj,jk)292 ! 342 ! ! now trend: 343 zu_frc(ji,jj) = zu_frc(ji,jj) + fse3u(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk) 344 zv_frc(ji,jj) = zv_frc(ji,jj) + fse3v(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk) 345 ! ! now bt transp: 346 un_b(ji,jj) = un_b(ji,jj) + fse3u(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk) 347 vn_b(ji,jj) = vn_b(ji,jj) + fse3v(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk) 348 ! ! before bt transp: 293 349 ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk) 294 350 vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk) … … 296 352 END DO 297 353 END DO 298 299 ! before inverse water column height at u- and v- points 354 ! 355 zu_frc(:,:) = zu_frc(:,:) * hur(:,:) 356 zv_frc(:,:) = zv_frc(:,:) * hvr(:,:) 357 ! 300 358 IF( lk_vvl ) THEN 301 zhu_b(:,:) = 0. 302 zhv_b(:,:) = 0. 303 DO jk = 1, jpk 304 zhu_b(:,:) = zhu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk) 305 zhv_b(:,:) = zhv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 306 END DO 307 zhu_b(:,:) = umask(:,:,1) / ( zhu_b(:,:) + 1. - umask(:,:,1) ) 308 zhv_b(:,:) = vmask(:,:,1) / ( zhv_b(:,:) + 1. - vmask(:,:,1) ) 309 ELSE 310 zhu_b(:,:) = hur(:,:) 311 zhv_b(:,:) = hvr(:,:) 312 ENDIF 313 359 ub_b(:,:) = ub_b(:,:) * zhur_b(:,:) 360 vb_b(:,:) = vb_b(:,:) * zhvr_b(:,:) 361 ELSE 362 ub_b(:,:) = ub_b(:,:) * hur(:,:) 363 vb_b(:,:) = vb_b(:,:) * hvr(:,:) 364 ENDIF 365 ! 314 366 ! !* baroclinic momentum trend (remove the vertical mean trend) 315 DO jk = 1, jpkm1 ! -------------------------- 367 DO jk = 1, jpkm1 ! ----------------------------------------------------------- 316 368 DO jj = 2, jpjm1 317 369 DO ji = fs_2, fs_jpim1 ! vector opt. 318 ua(ji,jj,jk) = ua(ji,jj,jk) - zu a(ji,jj) * hur(ji,jj)319 va(ji,jj,jk) = va(ji,jj,jk) - zv a(ji,jj) * hvr(ji,jj)370 ua(ji,jj,jk) = ua(ji,jj,jk) - zu_frc(ji,jj) * umask(ji,jj,jk) 371 va(ji,jj,jk) = va(ji,jj,jk) - zv_frc(ji,jj) * vmask(ji,jj,jk) 320 372 END DO 321 373 END DO 322 374 END DO 323 324 ! !* barotropic Coriolis trends * H (vorticity scheme dependent) 325 ! ! ---------------------------==== 326 zwx(:,:) = zun(:,:) * e2u(:,:) ! now transport 327 zwy(:,:) = zvn(:,:) * e1v(:,:) 375 ! !* barotropic Coriolis trends (vorticity scheme dependent) 376 ! ! -------------------------------------------------------- 377 zwx(:,:) = un_b(:,:) * e2u(:,:) ! now transport 378 zwy(:,:) = vn_b(:,:) * e1v(:,:) 328 379 ! 329 380 IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN ! energy conserving or mixed scheme … … 335 386 zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj) 336 387 ! energy conserving formulation for planetary vorticity term 337 z cu(ji,jj) = z1_4 * ( ff(ji ,jj-1) * zy1 + ff(ji,jj) * zy2 )338 z cv(ji,jj) =-z1_4 * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 )339 END DO 340 END DO 341 ! 342 ELSEIF ( ln_dynvor_ens ) THEN 388 zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 389 zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 390 END DO 391 END DO 392 ! 393 ELSEIF ( ln_dynvor_ens ) THEN ! enstrophy conserving scheme 343 394 DO jj = 2, jpjm1 344 395 DO ji = fs_2, fs_jpim1 ! vector opt. 345 zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj ) ) / e1u(ji,jj) 346 zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj) 347 zcu(ji,jj) = zy1 * ( ff(ji ,jj-1) + ff(ji,jj) ) 348 zcv(ji,jj) = zx1 * ( ff(ji-1,jj ) + ff(ji,jj) ) 349 END DO 350 END DO 351 ! 352 ELSEIF ( ln_dynvor_een ) THEN ! enstrophy and energy conserving scheme 396 zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 397 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj) 398 zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 399 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj) 400 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 401 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 402 END DO 403 END DO 404 ! 405 ELSEIF ( ln_dynvor_een ) THEN ! enstrophy and energy conserving scheme 353 406 DO jj = 2, jpjm1 354 407 DO ji = fs_2, fs_jpim1 ! vector opt. 355 zcu(ji,jj) = + z1_4 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 356 & + ftse(ji,jj ) * zwy(ji ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 357 zcv(ji,jj) = - z1_4 / e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji ,jj+1) & 358 & + ftnw(ji,jj ) * zwx(ji-1,jj ) + ftne(ji,jj ) * zwx(ji ,jj ) ) 359 END DO 360 END DO 361 ! 362 ENDIF 363 408 zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) & 409 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 410 & + ftse(ji,jj ) * zwy(ji ,jj-1) & 411 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 412 zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 413 & + ftse(ji,jj+1) * zwx(ji ,jj+1) & 414 & + ftnw(ji,jj ) * zwx(ji-1,jj ) & 415 & + ftne(ji,jj ) * zwx(ji ,jj ) ) 416 END DO 417 END DO 418 ! 419 ENDIF 420 ! 421 un_b (:,:) = un_b(:,:) * hur(:,:) ! Revert now transport to barotropic velocities 422 vn_b (:,:) = vn_b(:,:) * hvr(:,:) 364 423 ! !* Right-Hand-Side of the barotropic momentum equation 365 424 ! ! ---------------------------------------------------- 366 IF( lk_vvl ) THEN ! Variable volume : remove both Coriolis and Surface pressure gradient425 IF( lk_vvl ) THEN ! Variable volume : remove surface pressure gradient 367 426 DO jj = 2, jpjm1 368 427 DO ji = fs_2, fs_jpim1 ! vector opt. 369 zcu(ji,jj) = zcu(ji,jj) - grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn(ji+1,jj ) & 370 & - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) * hu(ji,jj) / e1u(ji,jj) 371 zcv(ji,jj) = zcv(ji,jj) - grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn(ji ,jj+1) & 372 & - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) * hv(ji,jj) / e2v(ji,jj) 373 END DO 374 END DO 375 ENDIF 376 377 DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend 428 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) / e1u(ji,jj) 429 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) / e2v(ji,jj) 430 END DO 431 END DO 432 ENDIF 433 434 DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend 378 435 DO ji = fs_2, fs_jpim1 379 zu a(ji,jj) = zua(ji,jj) - zcu(ji,jj)380 zv a(ji,jj) = zva(ji,jj) - zcv(ji,jj)436 zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * umask(ji,jj,1) 437 zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1) 381 438 END DO 382 END DO 383 384 385 ! ! Remove barotropic contribution of bottom friction 386 ! ! from the barotropic transport trend 387 zcoef = -1._wp * z1_2dt_b 388 389 IF(ln_bfrimp) THEN 390 ! ! Remove the bottom stress trend from 3-D sea surface level gradient 391 ! ! and Coriolis forcing in case of 3D semi-implicit bottom friction 392 DO jj = 2, jpjm1 393 DO ji = fs_2, fs_jpim1 394 ikbu = mbku(ji,jj) 395 ikbv = mbkv(ji,jj) 396 ua_btm = zcu(ji,jj) * z2dt_bf * hur(ji,jj) * umask (ji,jj,ikbu) 397 va_btm = zcv(ji,jj) * z2dt_bf * hvr(ji,jj) * vmask (ji,jj,ikbv) 398 399 zua(ji,jj) = zua(ji,jj) - bfrua(ji,jj) * ua_btm 400 zva(ji,jj) = zva(ji,jj) - bfrva(ji,jj) * va_btm 401 END DO 402 END DO 403 404 ELSE 405 406 # if defined key_vectopt_loop 407 DO jj = 1, 1 408 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling) 409 # else 410 DO jj = 2, jpjm1 411 DO ji = 2, jpim1 412 # endif 413 ! Apply stability criteria for bottom friction 414 !RBbug for vvl and external mode we may need to use varying fse3 415 !!gm Rq: the bottom e3 present the smallest variation, the use of e3u_0 is not a big approx. 416 zbfru(ji,jj) = MAX( bfrua(ji,jj) , fse3u(ji,jj,mbku(ji,jj)) * zcoef ) 417 zbfrv(ji,jj) = MAX( bfrva(ji,jj) , fse3v(ji,jj,mbkv(ji,jj)) * zcoef ) 418 END DO 419 END DO 420 421 DO jj = 2, jpjm1 422 DO ji = fs_2, fs_jpim1 ! vector opt. 423 zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) * zhu_b(ji,jj) 424 zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * zhv_b(ji,jj) 425 END DO 426 END DO 427 END IF ! end (ln_bfrimp) 428 429 430 ! !* d/dt(Ua), Ub, Vn (Vertical mean velocity) 431 ! ! -------------------------- 432 zua(:,:) = zua(:,:) * hur(:,:) 433 zva(:,:) = zva(:,:) * hvr(:,:) 434 ! 435 ub_b(:,:) = ub_b(:,:) * zhu_b(:,:) 436 vb_b(:,:) = vb_b(:,:) * zhv_b(:,:) 437 439 END DO 440 ! 441 ! ! Add bottom stress contribution from baroclinic velocities: 442 IF (ln_bt_fw) THEN 443 DO jj = 2, jpjm1 444 DO ji = fs_2, fs_jpim1 ! vector opt. 445 ikbu = mbku(ji,jj) 446 ikbv = mbkv(ji,jj) 447 zwx(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) ! NOW bottom baroclinic velocities 448 zwy(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj) 449 END DO 450 END DO 451 ELSE 452 DO jj = 2, jpjm1 453 DO ji = fs_2, fs_jpim1 ! vector opt. 454 ikbu = mbku(ji,jj) 455 ikbv = mbkv(ji,jj) 456 zwx(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) ! BEFORE bottom baroclinic velocities 457 zwy(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj) 458 END DO 459 END DO 460 ENDIF 461 ! 462 ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 463 zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:) 464 zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:) 465 ! 466 IF (ln_bt_fw) THEN ! Add wind forcing 467 zu_frc(:,:) = zu_frc(:,:) + zraur * utau(:,:) * hur(:,:) 468 zv_frc(:,:) = zv_frc(:,:) + zraur * vtau(:,:) * hvr(:,:) 469 ELSE 470 zu_frc(:,:) = zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * hur(:,:) 471 zv_frc(:,:) = zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * hvr(:,:) 472 ENDIF 473 ! 474 IF ( ln_apr_dyn ) THEN ! Add atm pressure forcing 475 IF (ln_bt_fw) THEN 476 DO jj = 2, jpjm1 477 DO ji = fs_2, fs_jpim1 ! vector opt. 478 zu_spg = grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) /e1u(ji,jj) 479 zv_spg = grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) /e2v(ji,jj) 480 zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg 481 zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg 482 END DO 483 END DO 484 ELSE 485 DO jj = 2, jpjm1 486 DO ji = fs_2, fs_jpim1 ! vector opt. 487 zu_spg = grav * z1_2 * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) & 488 & + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) /e1u(ji,jj) 489 zv_spg = grav * z1_2 * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) & 490 & + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) /e2v(ji,jj) 491 zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg 492 zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg 493 END DO 494 END DO 495 ENDIF 496 ENDIF 497 ! !* Right-Hand-Side of the barotropic ssh equation 498 ! ! ----------------------------------------------- 499 ! ! Surface net water flux and rivers 500 IF (ln_bt_fw) THEN 501 zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) ) 502 ELSE 503 zssh_frc(:,:) = zraur * z1_2 * (emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:)) 504 ENDIF 505 #if defined key_asminc 506 ! ! Include the IAU weighted SSH increment 507 IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN 508 zssh_frc(:,:) = zssh_frc(:,:) + ssh_iau(:,:) 509 ENDIF 510 #endif 511 ! 438 512 ! ----------------------------------------------------------------------- 439 ! Phase 2 : Integration of the barotropic equations with time splitting513 ! Phase 2 : Integration of the barotropic equations 440 514 ! ----------------------------------------------------------------------- 441 515 ! … … 456 530 ! Initialize depths: 457 531 IF ( lk_vvl.AND.(.NOT.ln_bt_fw) ) THEN 458 hu r_e (:,:) = zhu_b (:,:)459 hv r_e (:,:) = zhv_b (:,:)460 hu _e(:,:) = umask(:,:,1) / ( zhu_b(:,:) + 1. - umask(:,:,1))461 hv _e(:,:) = vmask(:,:,1) / ( zhv_b(:,:) + 1. - vmask(:,:,1))462 ELSE 463 hu_e 464 hv_e 465 hur_e 466 hvr_e 532 hu_e (:,:) = umask(:,:,1) / ( zhur_b(:,:) + 1._wp - umask(:,:,1) ) 533 hv_e (:,:) = vmask(:,:,1) / ( zhvr_b(:,:) + 1._wp - vmask(:,:,1) ) 534 hur_e (:,:) = zhur_b(:,:) 535 hvr_e (:,:) = zhvr_b(:,:) 536 ELSE 537 hu_e (:,:) = hu (:,:) 538 hv_e (:,:) = hv (:,:) 539 hur_e (:,:) = hur (:,:) 540 hvr_e (:,:) = hvr (:,:) 467 541 ENDIF 468 542 ! … … 479 553 zv_sum(:,:) = 0._wp 480 554 ! ! ==================== ! 481 482 zsshb_e(:,:) = sshn_b(:,:) ! sea surface height (before and now) 483 zub_e (:,:) = un_b (:,:) 484 zvb_e (:,:) = vn_b (:,:) 485 486 zu_sum (:,:) = un_b (:,:) ! summation 487 zv_sum (:,:) = vn_b (:,:) 488 zssh_sum(:,:) = sshn (:,:) 489 ua_b (:,:) = 0._wp ! After barotropic velocities (or transport if flux form) 490 va_b (:,:) = 0._wp 491 492 #if defined key_obc 493 ! set ssh corrections to 0 494 ! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop 495 IF( lp_obc_east ) sshfoe_b(:,:) = 0._wp 496 IF( lp_obc_west ) sshfow_b(:,:) = 0._wp 497 IF( lp_obc_south ) sshfos_b(:,:) = 0._wp 498 IF( lp_obc_north ) sshfon_b(:,:) = 0._wp 499 #endif 500 501 ! ! ==================== ! 502 DO jn = 1, icycle ! sub-time-step loop ! (from NOW to AFTER+1) 555 DO jn = 1, icycle ! sub-time-step loop ! 503 556 ! ! ==================== ! 504 z2dt_e = 2. * ( rdt / nn_baro )505 IF( jn == 1 ) z2dt_e = rdt / nn_baro506 507 557 ! !* Update the forcing (BDY and tides) 508 558 ! ! ------------------ 509 IF( lk_obc ) CALL obc_dta_bt( kt, jn ) 510 IF ( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset=(ioffset+1) ) 511 IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, koffset=ioffset ) 559 ! Update only tidal forcing at open boundaries 560 #if defined key_tide 561 IF ( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1) ) 562 IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, koffset=noffset ) 563 #endif 512 564 ! 513 565 ! Set extrapolation coefficients for predictor step: … … 523 575 524 576 ! Extrapolate barotropic velocities at step jit+0.5: 525 ua_e(:,:) = za1 * zun_e(:,:) + za2 * zub_e(:,:) + za3 * ubb_e(:,:)526 va_e(:,:) = za1 * zvn_e(:,:) + za2 * zvb_e(:,:) + za3 * vbb_e(:,:)577 ua_e(:,:) = za1 * zun_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:) 578 va_e(:,:) = za1 * zvn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:) 527 579 528 580 IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only) 529 581 ! ! ------------------ 530 582 ! Extrapolate Sea Level at step jit+0.5: 531 zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * zsshb_e(:,:) + za3 * sshbb_e(:,:)583 zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 532 584 ! 533 585 DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points … … 551 603 ! considering fluxes below: 552 604 ! 553 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)605 zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5 554 606 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) TIMESTEPS556 ! zwy(:,:) = e1v(:,:) * zvn_e(:,:) * hv_e(:,:)557 607 DO jj = 2, jpjm1 558 608 DO ji = fs_2, fs_jpim1 ! vector opt. … … 563 613 END DO 564 614 ! 565 #if defined key_obc 566 ! ! OBC : zhdiv must be zero behind the open boundary 567 !! mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column 568 IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1 ) = 0._wp ! east 569 IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1 ) = 0._wp ! west 570 IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0._wp ! north 571 IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1 ) = 0._wp ! south 572 #endif 573 #if defined key_bdy 574 zhdiv(:,:) = zhdiv(:,:) * bdytmask(:,:) ! BDY mask 575 #endif 576 ! 577 DO jj = 2, jpjm1 ! leap-frog on ssh_e 578 DO ji = fs_2, fs_jpim1 ! vector opt. 579 ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * ( zraur * ( emp(ji,jj)-rnf(ji,jj) ) + zhdiv(ji,jj) ) ) * tmask(ji,jj,1) 580 END DO 581 END DO 615 ! Sum over sub-time-steps to compute advective velocities 616 za2 = wgtbtp2(jn) 617 zu_sum (:,:) = zu_sum (:,:) + za2 * ua_e (:,:) * zhup2_e (:,:) 618 zv_sum (:,:) = zv_sum (:,:) + za2 * va_e (:,:) * zhvp2_e (:,:) 619 ! 620 ! Set next sea level: 621 ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1) 622 CALL lbc_lnk( ssha_e, 'T', 1._wp ) 623 582 624 #if defined key_bdy 583 625 ! Duplicate sea level across open boundaries (this is only cosmetic if lk_vvl=.false.) 584 626 IF (lk_bdy) CALL bdy_ssh( ssha_e ) 585 627 #endif 628 ! 629 ! Sea Surface Height at u-,v-points (vvl case only) 630 IF ( lk_vvl ) THEN 631 DO jj = 2, jpjm1 632 DO ji = 2, jpim1 ! NO Vector Opt. 633 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) / ( e1u(ji ,jj) * e2u(ji ,jj) ) & 634 & * ( e1t(ji ,jj) * e2t(ji ,jj) * ssha_e(ji ,jj) & 635 & + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha_e(ji+1,jj) ) 636 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) / ( e1v(ji,jj ) * e2v(ji,jj ) ) & 637 & * ( e1t(ji,jj ) * e2t(ji,jj ) * ssha_e(ji,jj ) & 638 & + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha_e(ji,jj+1) ) 639 END DO 640 END DO 641 CALL lbc_lnk( zsshu_a, 'U', 1._wp ) ; CALL lbc_lnk( zsshv_a, 'V', 1._wp ) 642 ENDIF 586 643 ! 587 644 ! Half-step back interpolation of SSH for surface pressure computation: … … 605 662 606 663 zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) & 607 & + za2 * zsshb_e(:,:) + za3 * sshbb_e(:,:)664 & + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 608 665 609 666 ! … … 635 692 DO jj = 2, jpjm1 636 693 DO ji = fs_2, fs_jpim1 ! vector opt. 637 ! surface pressure gradient638 IF( lk_vvl) THEN639 zu_spg = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) &640 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj)641 zv_spg = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) &642 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj)643 ELSE644 zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)645 zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)646 ENDIF647 ! add tidal astronomical forcing648 IF ( ln_tide_pot .AND. lk_tide ) THEN649 zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)650 zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)651 ENDIF652 ! energy conserving formulation for planetary vorticity term653 694 zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj) 654 695 zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj) 655 696 zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) / e2v(ji,jj) 656 697 zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj) 657 zu_cor = z1_4 * ( ff(ji ,jj-1) * zy1 + ff(ji,jj) * zy2 ) * hur_e(ji,jj) 658 zv_cor =-z1_4 * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj) 659 ! after velocities with implicit bottom friction 660 661 IF( ln_bfrimp ) THEN ! implicit bottom friction 662 ! A new method to implement the implicit bottom friction. 663 ! H. Liu 664 ! Sept 2011 665 ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) + & 666 & z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp ) & 667 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) 668 ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e * zua(ji,jj) ) * umask(ji,jj,1) 669 ! 670 va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) + & 671 & z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp ) & 672 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) 673 va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e * zva(ji,jj) ) * vmask(ji,jj,1) 674 ! 675 ELSE 676 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj))) * umask(ji,jj,1) & 677 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 678 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj))) * vmask(ji,jj,1) & 679 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 680 ENDIF 698 zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 699 zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 681 700 END DO 682 701 END DO … … 685 704 DO jj = 2, jpjm1 686 705 DO ji = fs_2, fs_jpim1 ! vector opt. 687 ! surface pressure gradient 688 IF( lk_vvl) THEN 689 zu_spg = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) & 690 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj) 691 zv_spg = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) & 692 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj) 693 ELSE 694 zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 695 zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 696 ENDIF 697 ! add tidal astronomical forcing 698 IF ( ln_tide_pot .AND. lk_tide ) THEN 699 zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 700 zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) 701 ENDIF 702 ! enstrophy conserving formulation for planetary vorticity term 703 zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj ) ) / e1u(ji,jj) 704 zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj) 705 zu_cor = zy1 * ( ff(ji ,jj-1) + ff(ji,jj) ) * hur_e(ji,jj) 706 zv_cor = zx1 * ( ff(ji-1,jj ) + ff(ji,jj) ) * hvr_e(ji,jj) 707 ! after velocities with implicit bottom friction 708 IF( ln_bfrimp ) THEN 709 ! A new method to implement the implicit bottom friction. 710 ! H. Liu 711 ! Sept 2011 712 ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) + & 713 & z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp ) & 714 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) 715 ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e * zua(ji,jj) ) * umask(ji,jj,1) 716 ! 717 va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) + & 718 & z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp ) & 719 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) 720 va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e * zva(ji,jj) ) * vmask(ji,jj,1) 721 ! 722 ELSE 723 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj))) * umask(ji,jj,1) & 724 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 725 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj))) * vmask(ji,jj,1) & 726 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 727 ENDIF 706 zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 707 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj) 708 zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 709 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj) 710 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 711 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 728 712 END DO 729 713 END DO … … 732 716 DO jj = 2, jpjm1 733 717 DO ji = fs_2, fs_jpim1 ! vector opt. 734 ! surface pressure gradient 735 IF( lk_vvl) THEN 736 zu_spg = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) & 737 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj) 738 zv_spg = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) & 739 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj) 740 ELSE 741 zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 742 zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 743 ENDIF 744 ! add tidal astronomical forcing 745 IF ( ln_tide_pot .AND. lk_tide ) THEN 746 zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 747 zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) 748 ENDIF 749 ! energy/enstrophy conserving formulation for planetary vorticity term 750 zu_cor = + z1_4 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 751 & + ftse(ji,jj ) * zwy(ji ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) * hur_e(ji,jj) 752 zv_cor = - z1_4 / e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji ,jj+1) & 753 & + ftnw(ji,jj ) * zwx(ji-1,jj ) + ftne(ji,jj ) * zwx(ji ,jj ) ) * hvr_e(ji,jj) 754 ! after velocities with implicit bottom friction 755 IF( ln_bfrimp ) THEN 756 ! A new method to implement the implicit bottom friction. 757 ! H. Liu 758 ! Sept 2011 759 ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) + & 760 & z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp ) & 761 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) 762 ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e * zua(ji,jj) ) * umask(ji,jj,1) 763 ! 764 va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) + & 765 & z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp ) & 766 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) 767 va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e * zva(ji,jj) ) * vmask(ji,jj,1) 768 ! 769 ELSE 770 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj))) * umask(ji,jj,1) & 771 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 772 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj))) * vmask(ji,jj,1) & 773 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 774 ENDIF 718 zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) & 719 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 720 & + ftse(ji,jj ) * zwy(ji ,jj-1) & 721 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 722 zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 723 & + ftse(ji,jj+1) * zwx(ji ,jj+1) & 724 & + ftnw(ji,jj ) * zwx(ji-1,jj ) & 725 & + ftne(ji,jj ) * zwx(ji ,jj ) ) 775 726 END DO 776 727 END DO 777 728 ! 778 729 ENDIF 779 ! !* domain lateral boundary 780 ! ! ----------------------- 781 782 ! OBC open boundaries 783 IF( lk_obc ) CALL obc_fla_ts ( ua_e, va_e, sshn_e, ssha_e ) 784 785 ! BDY open boundaries 786 #if defined key_bdy 787 pssh => sshn_e 730 ! 731 ! Add tidal astronomical forcing if defined 732 IF ( lk_tide.AND.ln_tide_pot ) THEN 733 DO jj = 2, jpjm1 734 DO ji = fs_2, fs_jpim1 ! vector opt. 735 zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 736 zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) 737 zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg 738 zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg 739 END DO 740 END DO 741 ENDIF 742 ! 743 ! Add bottom stresses: 744 zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:) 745 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:) 746 ! 747 ! Surface pressure trend: 748 DO jj = 2, jpjm1 749 DO ji = fs_2, fs_jpim1 ! vector opt. 750 ! Add surface pressure gradient 751 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 752 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 753 zwx(ji,jj) = zu_spg 754 zwy(ji,jj) = zv_spg 755 END DO 756 END DO 757 ! 758 ! Set next velocities: 759 IF( ln_dynadv_vec .OR. (.NOT. lk_vvl) ) THEN ! Vector form 760 DO jj = 2, jpjm1 761 DO ji = fs_2, fs_jpim1 ! vector opt. 762 ua_e(ji,jj) = ( zun_e(ji,jj) & 763 & + rdtbt * ( zwx(ji,jj) & 764 & + zu_trd(ji,jj) & 765 & + zu_frc(ji,jj) ) & 766 & ) * umask(ji,jj,1) 767 768 va_e(ji,jj) = ( zvn_e(ji,jj) & 769 & + rdtbt * ( zwy(ji,jj) & 770 & + zv_trd(ji,jj) & 771 & + zv_frc(ji,jj) ) & 772 & ) * vmask(ji,jj,1) 773 END DO 774 END DO 775 776 ELSE ! Flux form 777 DO jj = 2, jpjm1 778 DO ji = fs_2, fs_jpim1 ! vector opt. 779 780 zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1)) 781 zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1)) 782 783 ua_e(ji,jj) = ( hu_e(ji,jj) * zun_e(ji,jj) & 784 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 785 & + zhup2_e(ji,jj) * zu_trd(ji,jj) & 786 & + hu(ji,jj) * zu_frc(ji,jj) ) & 787 & ) * zhura 788 789 va_e(ji,jj) = ( hv_e(ji,jj) * zvn_e(ji,jj) & 790 & + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) & 791 & + zhvp2_e(ji,jj) * zv_trd(ji,jj) & 792 & + hv(ji,jj) * zv_frc(ji,jj) ) & 793 & ) * zhvra 794 END DO 795 END DO 796 ENDIF 797 ! 798 IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only) 799 ! ! ---------------------------------------------- 800 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 801 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 802 hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) ) 803 hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) ) 804 ! 805 ENDIF 806 ! !* domain lateral boundary 807 ! ! ----------------------- 808 ! 809 CALL lbc_lnk( ua_e , 'U', -1._wp ) ! local domain boundaries 810 CALL lbc_lnk( va_e , 'V', -1._wp ) 811 812 #if defined key_bdy 813 814 pssh => ssha_e 788 815 phur => hur_e 789 816 phvr => hvr_e 790 817 pu2d => ua_e 791 818 pv2d => va_e 792 793 IF( lk_bdy ) CALL bdy_dyn2d( kt ) 819 820 IF( lk_bdy ) CALL bdy_dyn2d( kt ) ! open boundaries 794 821 #endif 795 796 ! 797 CALL lbc_lnk( ua_e , 'U', -1. ) ! local domain boundaries 798 CALL lbc_lnk( va_e , 'V', -1. ) 799 CALL lbc_lnk( ssha_e, 'T', 1. ) 800 801 zu_sum (:,:) = zu_sum (:,:) + ua_e (:,:) ! Sum over sub-time-steps 802 zv_sum (:,:) = zv_sum (:,:) + va_e (:,:) 803 zssh_sum(:,:) = zssh_sum(:,:) + ssha_e(:,:) 804 805 ! !* Time filter and swap 806 ! ! -------------------- 807 sshbb_e(:,:) = zsshb_e(:,:) 808 ubb_e (:,:) = zub_e (:,:) 809 vbb_e (:,:) = zvb_e (:,:) 810 IF( jn == 1 ) THEN ! Swap only (1st Euler time step) 811 zsshb_e(:,:) = sshn_e(:,:) 812 zub_e (:,:) = zun_e (:,:) 813 zvb_e (:,:) = zvn_e (:,:) 814 sshn_e (:,:) = ssha_e(:,:) 815 zun_e (:,:) = ua_e (:,:) 816 zvn_e (:,:) = va_e (:,:) 817 ELSE ! Swap + Filter 818 zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:) 819 zub_e (:,:) = atfp * ( zub_e (:,:) + ua_e (:,:) ) + atfp1 * zun_e (:,:) 820 zvb_e (:,:) = atfp * ( zvb_e (:,:) + va_e (:,:) ) + atfp1 * zvn_e (:,:) 821 sshn_e (:,:) = ssha_e(:,:) 822 zun_e (:,:) = ua_e (:,:) 823 zvn_e (:,:) = va_e (:,:) 824 ENDIF 825 826 za1 = wgtbtp1(jn) 822 ! !* Swap 823 ! ! ---- 824 ubb_e (:,:) = ub_e (:,:) 825 ub_e (:,:) = zun_e (:,:) 826 zun_e (:,:) = ua_e (:,:) 827 ! 828 vbb_e (:,:) = vb_e (:,:) 829 vb_e (:,:) = zvn_e (:,:) 830 zvn_e (:,:) = va_e (:,:) 831 ! 832 sshbb_e(:,:) = sshb_e(:,:) 833 sshb_e (:,:) = sshn_e(:,:) 834 sshn_e (:,:) = ssha_e(:,:) 835 836 ! !* Sum over whole bt loop 837 ! ! ---------------------- 838 za1 = wgtbtp1(jn) 827 839 IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN ! Sum velocities 828 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) 829 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) 830 ELSE ! Sum transports 831 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) 832 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) 833 ! ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) 834 ! va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) 835 ENDIF 840 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) 841 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) 842 ELSE ! Sum transports 843 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) 844 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) 845 ENDIF 846 ! ! Sum sea level 836 847 ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:) 837 838 IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only)839 ! ! ------------------840 DO jj = 1, jpjm1 ! Sea Surface Height at u- & v-points841 DO ji = 1, fs_jpim1 ! Vector opt.842 zsshun_e(ji,jj) = 0.5_wp * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) &843 & * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn_e(ji ,jj) &844 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) )845 zsshvn_e(ji,jj) = 0.5_wp * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) &846 & * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn_e(ji,jj ) &847 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) )848 END DO849 END DO850 CALL lbc_lnk( zsshun_e, 'U', 1. ) ! lateral boundaries conditions851 CALL lbc_lnk( zsshvn_e, 'V', 1. )852 !853 hu_e (:,:) = hu_0(:,:) + zsshun_e(:,:) ! Ocean depth at U- and V-points854 hv_e (:,:) = hv_0(:,:) + zsshvn_e(:,:)855 hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) )856 hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) )857 !858 ENDIF859 848 ! ! ==================== ! 860 849 END DO ! end loop ! 861 850 ! ! ==================== ! 862 863 #if defined key_obc864 IF( lp_obc_east ) sshfoe_b(:,:) = zcoef * sshfoe_b(:,:) !!gm totally useless ?????865 IF( lp_obc_west ) sshfow_b(:,:) = zcoef * sshfow_b(:,:)866 IF( lp_obc_north ) sshfon_b(:,:) = zcoef * sshfon_b(:,:)867 IF( lp_obc_south ) sshfos_b(:,:) = zcoef * sshfos_b(:,:)868 #endif869 870 851 ! ----------------------------------------------------------------------------- 871 852 ! Phase 3. update the general trend with the barotropic trend 872 853 ! ----------------------------------------------------------------------------- 873 854 ! 874 ! !* Time average ==> after barotropic u, v, ssh 875 zcoef = 1._wp / ( icycle + 1 ) 876 zu_sum(:,:) = zcoef * zu_sum (:,:) 877 zv_sum(:,:) = zcoef * zv_sum (:,:) 878 ! 855 ! At this stage ssha holds a time averaged value 856 ! ! Sea Surface Height at u-,v- and f-points 857 IF( lk_vvl ) THEN ! (required only in key_vvl case) 858 DO jj = 1, jpjm1 859 DO ji = 1, jpim1 ! NO Vector Opt. 860 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) / ( e1u(ji ,jj) * e2u(ji ,jj) ) & 861 & * ( e1t(ji ,jj) * e2t(ji ,jj) * ssha(ji ,jj) & 862 & + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) ) 863 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) / ( e1v(ji,jj ) * e2v(ji,jj ) ) & 864 & * ( e1t(ji,jj ) * e2t(ji,jj ) * ssha(ji,jj ) & 865 & + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) ) 866 END DO 867 END DO 868 CALL lbc_lnk( zsshu_a, 'U', 1._wp ) ; CALL lbc_lnk( zsshv_a, 'V', 1._wp ) ! Boundary conditions 869 ENDIF 870 ! 879 871 ! Set advection velocity correction: 880 IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN 881 un_adv(:,:) = zu_sum(:,:) 882 vn_adv(:,:) = zv_sum(:,:) 883 ELSE 884 un_adv(:,:) = 0.5_wp * ( ub_b(:,:) + zu_sum(:,:))885 vn_adv(:,:) = 0.5_wp * ( vb_b(:,:) + zv_sum(:,:))872 IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN 873 un_adv(:,:) = zu_sum(:,:)*hur(:,:) 874 vn_adv(:,:) = zv_sum(:,:)*hvr(:,:) 875 ELSE 876 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zu_sum(:,:)) * hur(:,:) 877 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zv_sum(:,:)) * hvr(:,:) 886 878 END IF 887 ! !* update the general momentum trend 888 DO jk=1,jpkm1 889 ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) * z1_2dt_b 890 va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) * z1_2dt_b 891 END DO 892 un_b (:,:) = zu_sum(:,:) 893 vn_b (:,:) = zv_sum(:,:) 894 sshn_b(:,:) = zcoef * zssh_sum(:,:) 879 880 IF (ln_bt_fw) THEN ! Save integrated transport for next computation 881 ub2_b(:,:) = zu_sum(:,:) 882 vb2_b(:,:) = zv_sum(:,:) 883 ENDIF 884 ! 885 ! Update barotropic trend: 886 IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN 887 DO jk=1,jpkm1 888 ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b 889 va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b 890 END DO 891 ELSE 892 hu_e (:,:) = umask(:,:,1) / ( zhur_b(:,:) + 1._wp - umask(:,:,1) ) 893 hv_e (:,:) = vmask(:,:,1) / ( zhvr_b(:,:) + 1._wp - vmask(:,:,1) ) 894 DO jk=1,jpkm1 895 ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_e(:,:) ) * z1_2dt_b 896 va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_e(:,:) ) * z1_2dt_b 897 END DO 898 ! Save barotropic velocities not transport: 899 ua_b (:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - umask(:,:,1) ) 900 va_b (:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - vmask(:,:,1) ) 901 ENDIF 895 902 ! 896 903 ! !* write time-spliting arrays in the restart 897 IF( lrst_oce ) CALL ts_rst( kt, 'WRITE' ) 898 ! 899 CALL wrk_dealloc( jpi, jpj, zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv ) 900 CALL wrk_dealloc( jpi, jpj, zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e ) 901 CALL wrk_dealloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum ) 902 CALL wrk_dealloc( jpi, jpj, zhu_b, zhv_b ) 903 CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zsshp2_e ) 904 CALL wrk_dealloc( jpi, jpj, zhust_e, zhvst_e ) 904 IF(lrst_oce .AND.ln_bt_fw) CALL ts_rst( kt, 'WRITE' ) 905 ! 906 CALL wrk_dealloc( jpi, jpj, zsshp2_e, zhdiv ) 907 CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e ) 908 CALL wrk_dealloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc ) 909 CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e ) 910 CALL wrk_dealloc( jpi, jpj, zhur_b, zhvr_b ) 911 CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a ) 912 CALL wrk_dealloc( jpi, jpj, zht, zhf ) 905 913 ! 906 914 IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_ts') … … 994 1002 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 995 1003 ! 996 REAL(wp), POINTER, DIMENSION(:,:) :: zzhu_b, zzhv_b997 INTEGER :: ji, jk ! dummy loop indices998 1004 !!---------------------------------------------------------------------- 999 1005 ! 1000 1006 IF( TRIM(cdrw) == 'READ' ) THEN 1001 IF( iom_varid( numror, 'un_b', ldstop = .FALSE. ) > 0 ) THEN 1002 CALL iom_get( numror, jpdom_autoglo, 'un_b' , un_b (:,:) ) ! external velocity issued 1003 CALL iom_get( numror, jpdom_autoglo, 'vn_b' , vn_b (:,:) ) ! from barotropic loop 1004 ELSE 1005 un_b (:,:) = 0._wp 1006 vn_b (:,:) = 0._wp 1007 ! vertical sum 1008 IF( lk_vopt_loop ) THEN ! vector opt., forced unroll 1009 DO jk = 1, jpkm1 1010 DO ji = 1, jpij 1011 un_b(ji,1) = un_b(ji,1) + fse3u_n(ji,1,jk) * un(ji,1,jk) 1012 vn_b(ji,1) = vn_b(ji,1) + fse3v_n(ji,1,jk) * vn(ji,1,jk) 1013 END DO 1014 END DO 1015 ELSE ! No vector opt. 1016 DO jk = 1, jpkm1 1017 un_b(:,:) = un_b(:,:) + fse3u_n(:,:,jk) * un(:,:,jk) 1018 vn_b(:,:) = vn_b(:,:) + fse3v_n(:,:,jk) * vn(:,:,jk) 1019 END DO 1020 ENDIF 1021 un_b (:,:) = un_b(:,:) * hur(:,:) 1022 vn_b (:,:) = vn_b(:,:) * hvr(:,:) 1023 ENDIF 1024 1025 ! Vertically integrated velocity (before) 1026 IF (neuler/=0) THEN 1027 ub_b (:,:) = 0._wp 1028 vb_b (:,:) = 0._wp 1029 1030 ! vertical sum 1031 IF( lk_vopt_loop ) THEN ! vector opt., forced unroll 1032 DO jk = 1, jpkm1 1033 DO ji = 1, jpij 1034 ub_b(ji,1) = ub_b(ji,1) + fse3u_b(ji,1,jk) * ub(ji,1,jk) 1035 vb_b(ji,1) = vb_b(ji,1) + fse3v_b(ji,1,jk) * vb(ji,1,jk) 1036 END DO 1037 END DO 1038 ELSE ! No vector opt. 1039 DO jk = 1, jpkm1 1040 ub_b(:,:) = ub_b(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) 1041 vb_b(:,:) = vb_b(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) 1042 END DO 1043 ENDIF 1044 1045 IF( lk_vvl ) THEN 1046 CALL wrk_alloc( jpi, jpj, zzhu_b, zzhv_b ) 1047 ub_b (:,:) = 0. 1048 vb_b (:,:) = 0. 1049 zzhu_b(:,:) = 0. 1050 zzhv_b(:,:) = 0. 1051 ! vertical sum 1052 IF( lk_vopt_loop ) THEN ! vector opt., forced unroll 1053 DO jk = 1, jpkm1 1054 DO ji = 1, jpij 1055 ub_b (ji,1) = ub_b (ji,1) + fse3u_b(ji,1,jk) * ub (ji,1,jk) 1056 vb_b (ji,1) = vb_b (ji,1) + fse3v_b(ji,1,jk) * vb (ji,1,jk) 1057 zzhu_b(ji,1) = zzhu_b(ji,1) + fse3u_b(ji,1,jk) * umask(ji,1,jk) 1058 zzhv_b(ji,1) = zzhv_b(ji,1) + fse3v_b(ji,1,jk) * vmask(ji,1,jk) 1059 END DO 1060 END DO 1061 ELSE ! No vector opt. 1062 DO jk = 1, jpkm1 1063 ub_b (:,:) = ub_b (:,:) + fse3u_b(:,:,jk) * ub (:,:,jk) 1064 vb_b (:,:) = vb_b (:,:) + fse3v_b(:,:,jk) * vb (:,:,jk) 1065 zzhu_b(:,:) = zzhu_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk) 1066 zzhv_b(:,:) = zzhv_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 1067 END DO 1068 ENDIF 1069 ub_b(:,:) = ub_b(:,:) / ( zzhu_b(:,:) + 1. - umask(:,:,1) ) 1070 vb_b(:,:) = vb_b(:,:) / ( zzhv_b(:,:) + 1. - vmask(:,:,1) ) 1071 CALL wrk_dealloc( jpi, jpj, zzhu_b, zzhv_b ) 1072 ELSE 1073 ub_b (:,:) = 0.e0 1074 vb_b (:,:) = 0.e0 1075 ! vertical sum 1076 IF( lk_vopt_loop ) THEN ! vector opt., forced unroll 1077 DO jk = 1, jpkm1 1078 DO ji = 1, jpij 1079 ub_b(ji,1) = ub_b(ji,1) + fse3u_b(ji,1,jk) * ub(ji,1,jk) 1080 vb_b(ji,1) = vb_b(ji,1) + fse3v_b(ji,1,jk) * vb(ji,1,jk) 1081 END DO 1082 END DO 1083 ELSE ! No vector opt. 1084 DO jk = 1, jpkm1 1085 ub_b(:,:) = ub_b(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) 1086 vb_b(:,:) = vb_b(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) 1087 END DO 1088 ENDIF 1089 ub_b(:,:) = ub_b(:,:) * hur(:,:) 1090 vb_b(:,:) = vb_b(:,:) * hvr(:,:) 1091 ENDIF 1092 ELSE ! neuler==0 1093 ub_b (:,:) = un_b (:,:) 1094 vb_b (:,:) = vn_b (:,:) 1095 ENDIF 1096 1097 IF( iom_varid( numror, 'sshn_b', ldstop = .FALSE. ) > 0 ) THEN 1098 CALL iom_get( numror, jpdom_autoglo, 'sshn_b' , sshn_b (:,:) ) ! filtered ssh 1099 ELSE 1100 sshn_b(:,:) = sshb(:,:) ! if not in restart set previous time mean to current baroclinic before value 1101 ENDIF 1102 1007 CALL iom_get( numror, jpdom_autoglo, 'ub2_b' , ub2_b (:,:) ) 1008 CALL iom_get( numror, jpdom_autoglo, 'vb2_b' , vb2_b (:,:) ) 1103 1009 IF( .NOT.ln_bt_av .AND. iom_varid( numror, 'sshbb_e', ldstop = .FALSE. ) > 0) THEN 1104 1010 CALL iom_get( numror, jpdom_autoglo, 'sshbb_e' , sshbb_e(:,:) ) … … 1116 1022 vb_e = vb_b 1117 1023 ENDIF 1024 ! 1118 1025 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 1119 CALL iom_rstput( kt, nitrst, numrow, 'u n_b' , un_b (:,:) ) ! external velocity and ssh1120 CALL iom_rstput( kt, nitrst, numrow, 'v n_b' , vn_b (:,:) ) ! issued from barotropic loop1121 CALL iom_rstput( kt, nitrst, numrow, 'sshn_b' , sshn_b(:,:) ) !1026 CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:) ) 1027 CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:) ) 1028 ! 1122 1029 IF (.NOT.ln_bt_av) THEN 1123 1030 CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:) ) … … 1153 1060 CALL wrk_alloc( jpi, jpj, zcu ) 1154 1061 ! 1155 IF (lk_vvl) THEN 1156 DO jj = 1, jpj 1157 DO ji =1, jpi 1158 zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj)) 1159 zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj)) 1160 zcu(ji,jj) = sqrt(grav*ht_0(ji,jj)*(zxr2 + zyr2) ) 1161 END DO 1162 END DO 1163 ELSE 1164 DO jj = 1, jpj 1165 DO ji =1, jpi 1166 zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj)) 1167 zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj)) 1168 zcu(ji,jj) = sqrt( grav*ht_0(ji,jj)*(zxr2 + zyr2) ) 1169 END DO 1170 END DO 1171 ENDIF 1062 DO jj = 1, jpj 1063 DO ji =1, jpi 1064 zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj)) 1065 zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj)) 1066 zcu(ji,jj) = sqrt(grav*ht_0(ji,jj)*(zxr2 + zyr2) ) 1067 END DO 1068 END DO 1172 1069 1173 1070 zcmax = MAXVAL(zcu(:,:)) … … 1247 1144 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 1248 1145 WRITE(*,*) 'ts_rst : You should not have seen this print! error?', kt, cdrw 1249 END SUBROUTINE ts_rst 1146 END SUBROUTINE ts_rst 1250 1147 SUBROUTINE dyn_spg_ts_init( kt ) ! Empty routine 1251 1148 INTEGER , INTENT(in) :: kt ! ocean time-step
Note: See TracChangeset
for help on using the changeset viewer.