- Timestamp:
- 2013-11-20T17:28:04+01:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
r3764 r4292 17 17 !! 3.3 ! 2010-09 (D. Storkey, E.O'Dea) Bug fix for BDY module 18 18 !! 3.3 ! 2011-03 (P. Oddo) Bug fix for time-splitting+(BDY-OBC) and not VVL 19 !! 3.5 ! 2013-07 (J. Chanut) Compliant with time splitting changes 19 20 !!------------------------------------------------------------------------- 20 21 … … 42 43 USE wrk_nemo ! Memory Allocation 43 44 USE prtctl ! Print control 45 USE dynspg_ts ! Barotropic velocities 46 44 47 #if defined key_agrif 45 48 USE agrif_opa_interp … … 103 106 REAL(wp) :: zue3a, zue3n, zue3b, zuf, zec ! local scalars 104 107 REAL(wp) :: zve3a, zve3n, zve3b, zvf ! - - 108 REAL(wp), POINTER, DIMENSION(:,:) :: zua, zva 105 109 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3u_f, ze3v_f 106 110 !!---------------------------------------------------------------------- … … 109 113 ! 110 114 CALL wrk_alloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 115 IF ( lk_dynspg_ts ) CALL wrk_alloc( jpi,jpj, zua, zva ) 111 116 ! 112 117 IF( kt == nit000 ) THEN … … 127 132 ! 128 133 #else 134 135 # if defined key_dynspg_exp 129 136 ! Next velocity : Leap-frog time stepping 130 137 ! ------------- … … 147 154 END DO 148 155 ENDIF 149 156 # endif 157 158 # if defined key_dynspg_ts 159 ! Ensure below that barotropic velocities match time splitting estimate 160 ! Compute actual transport and replace it with ts estimate at "after" time step 161 zua(:,:) = 0._wp 162 zva(:,:) = 0._wp 163 IF (lk_vvl) THEN 164 DO jk = 1, jpkm1 165 zua(:,:) = zua(:,:) + fse3u_a(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 166 zva(:,:) = zva(:,:) + fse3v_a(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 167 END DO 168 DO jk = 1, jpkm1 169 ua(:,:,jk) = ( ua(:,:,jk) - zua(:,:) * hur_e(:,:) + ua_b(:,:) ) * umask(:,:,jk) 170 va(:,:,jk) = ( va(:,:,jk) - zva(:,:) * hvr_e(:,:) + va_b(:,:) ) * vmask(:,:,jk) 171 END DO 172 ELSE 173 DO jk = 1, jpkm1 174 zua(:,:) = zua(:,:) + fse3u(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 175 zva(:,:) = zva(:,:) + fse3v(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 176 END DO 177 DO jk = 1, jpkm1 178 ua(:,:,jk) = ( ua(:,:,jk) - zua(:,:) * hur(:,:) + ua_b(:,:) ) *umask(:,:,jk) 179 va(:,:,jk) = ( va(:,:,jk) - zva(:,:) * hvr(:,:) + va_b(:,:) ) *vmask(:,:,jk) 180 END DO 181 ENDIF 182 183 IF (lk_dynspg_ts.AND.(.NOT.ln_bt_fw)) THEN 184 ! Remove advective velocity from "now velocities" 185 ! prior to asselin filtering 186 ! In the forward case, this is done below after asselin filtering 187 DO jk = 1, jpkm1 188 un(:,:,jk) = ( un(:,:,jk) - un_adv(:,:) + un_b(:,:) )*umask(:,:,jk) 189 vn(:,:,jk) = ( vn(:,:,jk) - vn_adv(:,:) + vn_b(:,:) )*vmask(:,:,jk) 190 END DO 191 ENDIF 192 # endif 150 193 151 194 ! Update after velocity on domain lateral boundaries … … 194 237 vn(:,:,jk) = va(:,:,jk) 195 238 END DO 239 IF (lk_vvl) THEN 240 DO jk = 1, jpkm1 241 fse3t_b(:,:,jk) = fse3t_n(:,:,jk) 242 fse3u_b(:,:,jk) = fse3u_n(:,:,jk) 243 fse3v_b(:,:,jk) = fse3v_n(:,:,jk) 244 ENDDO 245 ENDIF 196 246 ELSE !* Leap-Frog : Asselin filter and swap 197 247 ! ! =============! … … 201 251 DO jj = 1, jpj 202 252 DO ji = 1, jpi 203 zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) )204 zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) )253 zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0_wp * un(ji,jj,jk) + ua(ji,jj,jk) ) 254 zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0_wp * vn(ji,jj,jk) + va(ji,jj,jk) ) 205 255 ! 206 256 ub(ji,jj,jk) = zuf ! ub <-- filtered velocity … … 214 264 ELSE ! Variable volume ! 215 265 ! ! ================! 266 ! Before scale factor at t-points 267 ! (used as a now filtered scale factor until the swap) 268 ! ---------------------------------------------------- 269 IF (lk_dynspg_ts.AND.ln_bt_fw) THEN 270 ! Remove asselin filtering on thicknesses if forward time splitting 271 fse3t_b(:,:,:) = fse3t_n(:,:,:) 272 ELSE 273 fse3t_b(:,:,:) = fse3t_n(:,:,:) + atfp * ( fse3t_b(:,:,:) - 2._wp * fse3t_n(:,:,:) + fse3t_a(:,:,:) ) 274 ! Add volume filter correction: compatibility with tracer advection scheme 275 ! => time filter + conservation correction (only at the first level) 276 fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 216 277 ! 217 DO jk = 1, jpkm1 ! Before scale factor at t-points 218 fse3t_b(:,:,jk) = fse3t_n(:,:,jk) & 219 & + atfp * ( fse3t_b(:,:,jk) + fse3t_a(:,:,jk) & 220 & - 2._wp * fse3t_n(:,:,jk) ) 221 END DO 222 zec = atfp * rdt / rau0 ! Add filter correction only at the 1st level of t-point scale factors 223 fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 278 ENDIF 224 279 ! 225 IF( ln_dynadv_vec ) THEN ! vector invariant form (no thickness weighted calulation) 226 ! 227 ! ! before scale factors at u- & v-pts (computed from fse3t_b) 228 CALL dom_vvl_2( kt, fse3u_b(:,:,:), fse3v_b(:,:,:) ) 229 ! 230 DO jk = 1, jpkm1 ! Leap-Frog - Asselin filter and swap: applied on velocity 231 DO jj = 1, jpj ! -------- 280 IF( ln_dynadv_vec ) THEN 281 ! Before scale factor at (u/v)-points 282 ! ----------------------------------- 283 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' ) 284 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' ) 285 ! Leap-Frog - Asselin filter and swap: applied on velocity 286 ! ----------------------------------- 287 DO jk = 1, jpkm1 288 DO jj = 1, jpj 232 289 DO ji = 1, jpi 233 zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2. e0* un(ji,jj,jk) + ua(ji,jj,jk) )234 zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2. e0* vn(ji,jj,jk) + va(ji,jj,jk) )290 zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2._wp * un(ji,jj,jk) + ua(ji,jj,jk) ) 291 zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2._wp * vn(ji,jj,jk) + va(ji,jj,jk) ) 235 292 ! 236 293 ub(ji,jj,jk) = zuf ! ub <-- filtered velocity … … 242 299 END DO 243 300 ! 244 ELSE ! flux form (thickness weighted calulation) 245 ! 246 CALL dom_vvl_2( kt, ze3u_f, ze3v_f ) ! before scale factors at u- & v-pts (computed from fse3t_b) 247 ! 248 DO jk = 1, jpkm1 ! Leap-Frog - Asselin filter and swap: 249 DO jj = 1, jpj ! applied on thickness weighted velocity 301 ELSE 302 ! Temporary filtered scale factor at (u/v)-points (will become before scale factor) 303 !------------------------------------------------ 304 CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3u_f, 'U' ) 305 CALL dom_vvl_interpol( fse3t_b(:,:,:), ze3v_f, 'V' ) 306 ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity 307 ! ----------------------------------- =========================== 308 DO jk = 1, jpkm1 309 DO jj = 1, jpj 250 310 DO ji = 1, jpi ! --------------------------- 251 311 zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk) … … 272 332 ENDIF 273 333 ! 274 ENDIF 334 IF (lk_dynspg_ts.AND.ln_bt_fw) THEN 335 ! Remove asselin filtering of barotropic velocities if forward time splitting 336 ! note that we replace barotropic velocities by advective velocities 337 zua(:,:) = 0._wp 338 zva(:,:) = 0._wp 339 IF (lk_vvl) THEN 340 DO jk = 1, jpkm1 341 zua(:,:) = zua(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) * umask(:,:,jk) 342 zva(:,:) = zva(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk) 343 END DO 344 ELSE 345 DO jk = 1, jpkm1 346 zua(:,:) = zua(:,:) + fse3u(:,:,jk) * ub(:,:,jk) * umask(:,:,jk) 347 zva(:,:) = zva(:,:) + fse3v(:,:,jk) * vb(:,:,jk) * vmask(:,:,jk) 348 END DO 349 ENDIF 350 DO jk = 1, jpkm1 351 ub(:,:,jk) = ub(:,:,jk) - (zua(:,:) * hur(:,:) - un_b(:,:)) * umask(:,:,jk) 352 vb(:,:,jk) = vb(:,:,jk) - (zva(:,:) * hvr(:,:) - vn_b(:,:)) * vmask(:,:,jk) 353 END DO 354 ENDIF 355 ! 356 ENDIF ! neuler =/0 275 357 276 358 IF(ln_ctl) CALL prt_ctl( tab3d_1=un, clinfo1=' nxt - Un: ', mask1=umask, & … … 278 360 ! 279 361 CALL wrk_dealloc( jpi,jpj,jpk, ze3u_f, ze3v_f ) 362 IF ( lk_dynspg_ts ) CALL wrk_dealloc( jpi,jpj, zua, zva ) 280 363 ! 281 364 IF( nn_timing == 1 ) CALL timing_stop('dyn_nxt')
Note: See TracChangeset
for help on using the changeset viewer.
