- Timestamp:
- 2015-12-16T10:25:22+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
r5930 r6060 2 2 !!====================================================================== 3 3 !! *** MODULE dynzdf_imp *** 4 !! Ocean dynamics: vertical component(s) of the momentum mixing trend 4 !! Ocean dynamics: vertical component(s) of the momentum mixing trend, implicit scheme 5 5 !!====================================================================== 6 6 !! History : OPA ! 1990-10 (B. Blanke) Original code … … 12 12 13 13 !!---------------------------------------------------------------------- 14 !! dyn_zdf_imp : update the momentum trend with the vertical diffusion using a implicit time-stepping 15 !!---------------------------------------------------------------------- 16 USE oce ! ocean dynamics and tracers 17 USE dom_oce ! ocean space and time domain 18 USE domvvl ! variable volume 19 USE sbc_oce ! surface boundary condition: ocean 20 USE zdf_oce ! ocean vertical physics 21 USE phycst ! physical constants 22 USE in_out_manager ! I/O manager 23 USE lib_mpp ! MPP library 24 USE zdfbfr ! Bottom friction setup 25 USE wrk_nemo ! Memory Allocation 26 USE timing ! Timing 27 USE dynadv, ONLY: ln_dynadv_vec ! Momentum advection form 14 !! dyn_zdf_imp : compute the vertical diffusion using a implicit scheme 15 !! together with the Leap-Frog time integration. 16 !!---------------------------------------------------------------------- 17 USE oce ! ocean dynamics and tracers 18 USE phycst ! physical constants 19 USE dom_oce ! ocean space and time domain 20 USE domvvl ! variable volume 21 USE sbc_oce ! surface boundary condition: ocean 22 USE dynadv , ONLY: ln_dynadv_vec ! Momentum advection form 23 USE zdf_oce ! ocean vertical physics 24 USE zdfbfr ! Bottom friction setup 25 ! 26 USE in_out_manager ! I/O manager 27 USE lib_mpp ! MPP library 28 USE wrk_nemo ! Memory Allocation 29 USE timing ! Timing 28 30 29 31 IMPLICIT NONE … … 32 34 PUBLIC dyn_zdf_imp ! called by step.F90 33 35 34 REAL(wp) :: r_vvl ! variable volume indicator, =1 if lk_vvl=T, =0otherwise36 REAL(wp) :: r_vvl ! non-linear free surface indicator: =0 if ln_linssh=T, =1 otherwise 35 37 36 38 !! * Substitutions 37 # include "domzgr_substitute.h90"38 39 # include "vectopt_loop_substitute.h90" 39 40 !!---------------------------------------------------------------------- … … 49 50 !! 50 51 !! ** Purpose : Compute the trend due to the vert. momentum diffusion 51 !! and the surface forcing, and add it to the general trend of52 !! the momentum equations.52 !! together with the Leap-Frog time stepping using an 53 !! implicit scheme. 53 54 !! 54 !! ** Method : The vertical momentum mixing trend is given by : 55 !! dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) ) 56 !! backward time stepping 57 !! Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2) 58 !! Bottom boundary conditions : bottom stress (cf zdfbfr.F) 59 !! Add this trend to the general trend ua : 60 !! ua = ua + dz( avmu dz(u) ) 55 !! ** Method : - Leap-Frog time stepping on all trends but the vertical mixing 56 !! ua = ub + 2*dt * ua vector form or linear free surf. 57 !! ua = ( e3u_b*ub + 2*dt * e3u_n*ua ) / e3u_a otherwise 58 !! - update the after velocity with the implicit vertical mixing. 59 !! This requires to solver the following system: 60 !! ua = ua + 1/e3u_a dk+1[ avmu / e3uw_a dk[ua] ] 61 !! with the following surface/top/bottom boundary condition: 62 !! surface: wind stress input (averaged over kt-1/2 & kt+1/2) 63 !! top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfbfr.F) 61 64 !! 62 !! ** Action : - Update (ua,va) arrays with the after vertical diffusive mixing trend.65 !! ** Action : (ua,va) after velocity 63 66 !!--------------------------------------------------------------------- 64 67 INTEGER , INTENT(in) :: kt ! ocean time-step index 65 68 REAL(wp), INTENT(in) :: p2dt ! vertical profile of tracer time-step 66 ! !67 INTEGER :: ji, jj, jk ! dummy loop indices68 INTEGER :: ikbu, ikbv ! local integers69 REAL(wp) :: z 1_p2dt, zcoef, zzwi, zzws, zrhs! local scalars70 REAL(wp) :: z e3ua, ze3va69 ! 70 INTEGER :: ji, jj, jk ! dummy loop indices 71 INTEGER :: ikbu, ikbv ! local integers 72 REAL(wp) :: zzwi, ze3ua ! local scalars 73 REAL(wp) :: zzws, ze3va ! - - 71 74 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwd, zws 72 75 !!---------------------------------------------------------------------- … … 81 84 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 82 85 ! 83 I F( lk_vvl ) THEN ; r_vvl = 1._wp ! Variable volume indicator84 ELSE ; r_vvl = 0._wp86 If( ln_linssh ) THEN ; r_vvl = 0._wp ! non-linear free surface indicator 87 ELSE ; r_vvl = 1._wp 85 88 ENDIF 86 89 ENDIF 87 88 ! 1. Time step dynamics 89 ! --------------------- 90 91 IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN ! applied on velocity 90 ! 91 ! !== Time step dynamics ==! 92 ! 93 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity 92 94 DO jk = 1, jpkm1 93 95 ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk) 94 96 va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk) 95 97 END DO 96 ELSE 98 ELSE ! applied on thickness weighted velocity 97 99 DO jk = 1, jpkm1 98 ua(:,:,jk) = ( ub(:,:,jk) * fse3u_b(:,:,jk) & 99 & + p2dt * ua(:,:,jk) * fse3u_n(:,:,jk) ) & 100 & / fse3u_a(:,:,jk) * umask(:,:,jk) 101 va(:,:,jk) = ( vb(:,:,jk) * fse3v_b(:,:,jk) & 102 & + p2dt * va(:,:,jk) * fse3v_n(:,:,jk) ) & 103 & / fse3v_a(:,:,jk) * vmask(:,:,jk) 104 END DO 105 ENDIF 106 107 ! 2. Apply semi-implicit bottom friction 108 ! -------------------------------------- 100 ua(:,:,jk) = ( e3u_b(:,:,jk) * ub(:,:,jk) & 101 & + p2dt * e3u_n(:,:,jk) * ua(:,:,jk) ) / e3u_a(:,:,jk) * umask(:,:,jk) 102 va(:,:,jk) = ( e3v_b(:,:,jk) * vb(:,:,jk) & 103 & + p2dt * e3v_n(:,:,jk) * va(:,:,jk) ) / e3v_a(:,:,jk) * vmask(:,:,jk) 104 END DO 105 ENDIF 106 ! 107 ! !== Apply semi-implicit bottom friction ==! 108 ! 109 109 ! Only needed for semi-implicit bottom friction setup. The explicit 110 110 ! bottom friction has been included in "u(v)a" which act as the R.H.S … … 116 116 ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points 117 117 ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 118 avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1)119 avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1)118 avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * e3uw_n(ji,jj,ikbu+1) 119 avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * e3vw_n(ji,jj,ikbv+1) 120 120 END DO 121 121 END DO … … 125 125 ikbu = miku(ji,jj) ! ocean top level at u- and v-points 126 126 ikbv = mikv(ji,jj) ! (first wet ocean u- and v-points) 127 IF (ikbu .GE. 2) avmu(ji,jj,ikbu) = -tfrua(ji,jj) * fse3uw(ji,jj,ikbu)128 IF (ikbv .GE. 2) avmv(ji,jj,ikbv) = -tfrva(ji,jj) * fse3vw(ji,jj,ikbv)127 IF( ikbu >= 2 ) avmu(ji,jj,ikbu) = -tfrua(ji,jj) * e3uw_n(ji,jj,ikbu) 128 IF( ikbv >= 2 ) avmv(ji,jj,ikbv) = -tfrva(ji,jj) * e3vw_n(ji,jj,ikbv) 129 129 END DO 130 130 END DO 131 131 END IF 132 132 ENDIF 133 133 ! 134 134 ! With split-explicit free surface, barotropic stress is treated explicitly 135 135 ! Update velocities at the bottom. 136 136 ! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does 137 137 ! not lead to the effective stress seen over the whole barotropic loop. 138 IF ( ln_bfrimp .AND.ln_dynspg_ts ) THEN 139 ! remove barotropic velocities: 140 DO jk = 1, jpkm1 141 ua(:,:,jk) = (ua(:,:,jk) - ua_b(:,:)) * umask(:,:,jk) 142 va(:,:,jk) = (va(:,:,jk) - va_b(:,:)) * vmask(:,:,jk) 143 END DO 144 ! Add bottom/top stress due to barotropic component only: 145 DO jj = 2, jpjm1 138 ! G. Madec : in linear free surface, e3u_a = e3u_n = e3u_0, so systematic use of e3u_a 139 IF( ln_bfrimp .AND. ln_dynspg_ts ) THEN 140 DO jk = 1, jpkm1 ! remove barotropic velocities 141 ua(:,:,jk) = ( ua(:,:,jk) - ua_b(:,:) ) * umask(:,:,jk) 142 va(:,:,jk) = ( va(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk) 143 END DO 144 DO jj = 2, jpjm1 ! Add bottom/top stress due to barotropic component only 146 145 DO ji = fs_2, fs_jpim1 ! vector opt. 147 146 ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points 148 147 ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 149 ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl * fse3u_a(ji,jj,ikbu)150 ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl * fse3v_a(ji,jj,ikbv)148 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,ikbu) + r_vvl * e3u_a(ji,jj,ikbu) 149 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikbv) + r_vvl * e3v_a(ji,jj,ikbv) 151 150 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * bfrua(ji,jj) * ua_b(ji,jj) / ze3ua 152 151 va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * bfrva(ji,jj) * va_b(ji,jj) / ze3va 153 152 END DO 154 153 END DO 155 IF ( ln_isfcav ) THEN154 IF( ln_isfcav ) THEN ! Ocean cavities (ISF) 156 155 DO jj = 2, jpjm1 157 156 DO ji = fs_2, fs_jpim1 ! vector opt. 158 157 ikbu = miku(ji,jj) ! top ocean level at u- and v-points 159 158 ikbv = mikv(ji,jj) ! (first wet ocean u- and v-points) 160 ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,ikbu) + r_vvl * fse3u_a(ji,jj,ikbu)161 ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj,ikbv) + r_vvl * fse3v_a(ji,jj,ikbv)159 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,ikbu) + r_vvl * e3u_a(ji,jj,ikbu) 160 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikbv) + r_vvl * e3v_a(ji,jj,ikbv) 162 161 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * tfrua(ji,jj) * ua_b(ji,jj) / ze3ua 163 162 va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * tfrva(ji,jj) * va_b(ji,jj) / ze3va … … 166 165 END IF 167 166 ENDIF 168 169 ! 3. Vertical diffusion on u170 ! ---------------------------167 ! 168 ! !== Vertical diffusion on u ==! 169 ! 171 170 ! Matrix and second member construction 172 171 ! bottom boundary condition: both zwi and zws must be masked as avmu can take … … 176 175 DO jj = 2, jpjm1 177 176 DO ji = fs_2, fs_jpim1 ! vector opt. 178 ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,jk) + r_vvl * fse3u_a(ji,jj,jk) ! after scale factor at T-point 179 zcoef = - p2dt / ze3ua 180 zzwi = zcoef * avmu (ji,jj,jk ) / fse3uw(ji,jj,jk ) 181 zwi(ji,jj,jk) = zzwi * wumask(ji,jj,jk ) 182 zzws = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 183 zws(ji,jj,jk) = zzws * wumask(ji,jj,jk+1) 177 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at T-point 178 zzwi = - p2dt * avmu(ji,jj,jk ) / ( ze3ua * e3uw_n(ji,jj,jk ) ) 179 zzws = - p2dt * avmu(ji,jj,jk+1) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) 180 zwi(ji,jj,jk) = zzwi * wumask(ji,jj,jk ) 181 zws(ji,jj,jk) = zzws * wumask(ji,jj,jk+1) 184 182 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 185 183 END DO … … 216 214 END DO 217 215 ! 218 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 219 DO ji = fs_2, fs_jpim1 ! vector opt. 220 ze3ua = ( 1._wp - r_vvl ) * fse3u_n(ji,jj,1) + r_vvl * fse3u_a(ji,jj,1)216 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==! 217 DO ji = fs_2, fs_jpim1 ! vector opt. 218 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 221 219 ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 222 220 & / ( ze3ua * rau0 ) * umask(ji,jj,1) … … 226 224 DO jj = 2, jpjm1 227 225 DO ji = fs_2, fs_jpim1 228 zrhs = ua(ji,jj,jk) ! zrhs=right hand side 229 ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 230 END DO 231 END DO 232 END DO 233 ! 234 DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk == 226 ua(ji,jj,jk) = ua(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 227 END DO 228 END DO 229 END DO 230 ! 231 DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk ==! 235 232 DO ji = fs_2, fs_jpim1 ! vector opt. 236 233 ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) … … 244 241 END DO 245 242 END DO 246 247 ! 4. Vertical diffusion on v248 ! ---------------------------243 ! 244 ! !== Vertical diffusion on v ==! 245 ! 249 246 ! Matrix and second member construction 250 247 ! bottom boundary condition: both zwi and zws must be masked as avmv can take … … 254 251 DO jj = 2, jpjm1 255 252 DO ji = fs_2, fs_jpim1 ! vector opt. 256 ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj,jk) + r_vvl * fse3v_a(ji,jj,jk) ! after scale factor at T-point 257 zcoef = - p2dt / ze3va 258 zzwi = zcoef * avmv (ji,jj,jk ) / fse3vw(ji,jj,jk ) 259 zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk) 260 zzws = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1) 261 zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) 253 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at T-point 254 zzwi = - p2dt * avmv (ji,jj,jk ) / ( ze3va * e3vw_n(ji,jj,jk ) ) 255 zzws = - p2dt * avmv (ji,jj,jk+1) / ( ze3va * e3vw_n(ji,jj,jk+1) ) 256 zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk ) 257 zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) 262 258 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 263 259 END DO … … 294 290 END DO 295 291 ! 296 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 292 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==! 297 293 DO ji = fs_2, fs_jpim1 ! vector opt. 298 ze3va = ( 1._wp - r_vvl ) * fse3v_n(ji,jj,1) + r_vvl * fse3v_a(ji,jj,1)294 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 299 295 va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 300 296 & / ( ze3va * rau0 ) … … 304 300 DO jj = 2, jpjm1 305 301 DO ji = fs_2, fs_jpim1 ! vector opt. 306 zrhs = va(ji,jj,jk) ! zrhs=right hand side 307 va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 308 END DO 309 END DO 310 END DO 311 ! 312 DO jj = 2, jpjm1 !== third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk == 302 va(ji,jj,jk) = va(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 303 END DO 304 END DO 305 END DO 306 ! 307 DO jj = 2, jpjm1 !== third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk ==! 313 308 DO ji = fs_2, fs_jpim1 ! vector opt. 314 309 va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) … … 322 317 END DO 323 318 END DO 324 319 325 320 ! J. Chanut: Lines below are useless ? 326 !! restore avmu(v)=0. at bottom (and top if ln_isfcav=T interfaces) 321 !! restore bottom layer avmu(v) 322 !!gm I almost sure it is !!!! 327 323 IF( ln_bfrimp ) THEN 328 324 DO jj = 2, jpjm1 … … 330 326 ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points 331 327 ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 332 avmu(ji,jj,ikbu+1) = 0. e0333 avmv(ji,jj,ikbv+1) = 0. e0328 avmu(ji,jj,ikbu+1) = 0._wp 329 avmv(ji,jj,ikbv+1) = 0._wp 334 330 END DO 335 331 END DO … … 339 335 ikbu = miku(ji,jj) ! ocean top level at u- and v-points 340 336 ikbv = mikv(ji,jj) ! (first wet ocean u- and v-points) 341 IF (ikbu > 1) avmu(ji,jj,ikbu) = 0.e0342 IF (ikbv > 1) avmv(ji,jj,ikbv) = 0.e0337 IF( ikbu > 1 ) avmu(ji,jj,ikbu) = 0._wp 338 IF( ikbv > 1 ) avmv(ji,jj,ikbv) = 0._wp 343 339 END DO 344 340 END DO 345 END 346 ENDIF 347 ! 348 CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwd, zws)349 ! 350 IF( nn_timing == 1 ) CALL timing_stop('dyn_zdf_imp')341 ENDIF 342 ENDIF 343 ! 344 CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwd, zws) 345 ! 346 IF( nn_timing == 1 ) CALL timing_stop('dyn_zdf_imp') 351 347 ! 352 348 END SUBROUTINE dyn_zdf_imp
Note: See TracChangeset
for help on using the changeset viewer.
