- Timestamp:
- 2017-06-10T09:31:34+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r7881_HPC09_ZDF/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90
r7990 r8160 6 6 !! History : 1.0 ! 2005-11 (G. Madec) Original code 7 7 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase 8 !! 4.0 ! 2017-06 (G. Madec) remove the explicit time-stepping option 8 9 !!---------------------------------------------------------------------- 9 10 10 11 !!---------------------------------------------------------------------- 11 12 !! dyn_zdf : Update the momentum trend with the vertical diffusion 12 !! dyn_zdf_init : initializations of the vertical diffusion scheme13 13 !!---------------------------------------------------------------------- 14 14 USE oce ! ocean dynamics and tracers variables 15 USE phycst ! physical constants 15 16 USE dom_oce ! ocean space and time domain variables 17 USE sbc_oce ! surface boundary condition: ocean 16 18 USE zdf_oce ! ocean vertical physics variables 17 USE dynzdf_exp ! vertical diffusion: explicit (dyn_zdf_exp routine)18 USE dyn zdf_imp ! vertical diffusion: implicit (dyn_zdf_imp routine)19 USE zdfdrg ! vertical physics: top/bottom drag coef. 20 USE dynadv , ONLY: ln_dynadv_vec ! dynamics: advection form 19 21 USE ldfdyn ! lateral diffusion: eddy viscosity coef. 20 22 USE trd_oce ! trends: ocean variables … … 30 32 PRIVATE 31 33 32 PUBLIC dyn_zdf ! routine called by step.F90 33 PUBLIC dyn_zdf_init ! routine called by opa.F90 34 35 INTEGER :: nzdf = 0 ! type vertical diffusion algorithm used, defined from ln_zdf... namlist logicals 34 PUBLIC dyn_zdf ! routine called by step.F90 35 36 REAL(wp) :: r_vvl ! non-linear free surface indicator: =0 if ln_linssh=T, =1 otherwise 36 37 37 38 !! * Substitutions 38 39 # include "vectopt_loop_substitute.h90" 39 40 !!---------------------------------------------------------------------- 40 !! NEMO/OPA 3.3 , NEMO Consortium (2010)41 !! NEMO/OPA 4.0 , NEMO Consortium (2017) 41 42 !! $Id$ 42 43 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 57 58 IF( nn_timing == 1 ) CALL timing_start('dyn_zdf') 58 59 ! 59 ! !set time step60 ! !* set time step 60 61 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt = rdt ! = rdt (restart with Euler time stepping) 61 62 ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt = 2. * rdt ! = 2 rdt (leapfrog) 62 63 ENDIF 63 64 64 IF( l_trddyn ) THEN !temporary save of ta and sa trends65 IF( l_trddyn ) THEN !* temporary save of ta and sa trends 65 66 CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 66 67 ztrdu(:,:,:) = ua(:,:,:) 67 68 ztrdv(:,:,:) = va(:,:,:) 68 69 ENDIF 69 70 SELECT CASE ( nzdf ) ! compute lateral mixing trend and add it to the general trend 71 ! 72 CASE ( 0 ) ; CALL dyn_zdf_exp( kt, r2dt ) ! explicit scheme 73 CASE ( 1 ) ; CALL dyn_zdf_imp( kt, r2dt ) ! implicit scheme 74 ! 75 END SELECT 76 70 ! !* compute lateral mixing trend and add it to the general trend 71 ! 72 CALL dyn_zdf_imp( kt, r2dt ) 73 ! 77 74 IF( l_trddyn ) THEN ! save the vertical diffusive trends for further diagnostics 78 75 ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / r2dt - ztrdu(:,:,:) … … 90 87 91 88 92 SUBROUTINE dyn_zdf_i nit89 SUBROUTINE dyn_zdf_imp( kt, p2dt ) 93 90 !!---------------------------------------------------------------------- 94 !! *** ROUTINE dyn_zdf_init *** 91 !! *** ROUTINE dyn_zdf_imp *** 92 !! 93 !! ** Purpose : Compute the trend due to the vert. momentum diffusion 94 !! together with the Leap-Frog time stepping using an 95 !! implicit scheme. 95 96 !! 96 !! ** Purpose : initialization of the vertical diffusion scheme 97 !! ** Method : - Leap-Frog time stepping on all trends but the vertical mixing 98 !! ua = ub + 2*dt * ua vector form or linear free surf. 99 !! ua = ( e3u_b*ub + 2*dt * e3u_n*ua ) / e3u_a otherwise 100 !! - update the after velocity with the implicit vertical mixing. 101 !! This requires to solver the following system: 102 !! ua = ua + 1/e3u_a dk+1[ avmu / e3uw_a dk[ua] ] 103 !! with the following surface/top/bottom boundary condition: 104 !! surface: wind stress input (averaged over kt-1/2 & kt+1/2) 105 !! top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfdrg.F90) 97 106 !! 98 !! ** Method : implicit (euler backward) scheme (default) 99 !! explicit (time-splitting) scheme if ln_zdfexp=T 107 !! ** Action : (ua,va) after velocity 108 !!--------------------------------------------------------------------- 109 INTEGER , INTENT(in) :: kt ! ocean time-step index 110 REAL(wp), INTENT(in) :: p2dt ! vertical profile of tracer time-step 111 ! 112 INTEGER :: ji, jj, jk ! dummy loop indices 113 INTEGER :: ikbu, ikbv ! local integers 114 REAL(wp) :: zzwi, ze3ua ! local scalars 115 REAL(wp) :: zzws, ze3va ! - - 116 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi, zwd, zws 100 117 !!---------------------------------------------------------------------- 101 USE zdftke 102 USE zdfgls 103 !!---------------------------------------------------------------------- 104 ! 105 ! Choice from ln_zdfexp (namzdf namelist variable read in zdfphy module) 106 IF( ln_zdfexp ) THEN ; nzdf = 0 ! use explicit scheme 107 ELSE ; nzdf = 1 ! use implicit scheme 108 ENDIF 109 ! 110 ! Force implicit schemes 111 IF( ln_zdftke .OR. ln_zdfgls ) nzdf = 1 ! TKE or GLS physics 112 IF( ln_dynldf_iso ) nzdf = 1 ! iso-neutral lateral physics 113 IF( ln_dynldf_hor .AND. ln_sco ) nzdf = 1 ! horizontal lateral physics in s-coordinate 114 ! 115 IF(lwp) THEN ! Print the choice 116 WRITE(numout,*) 117 WRITE(numout,*) 'dyn_zdf_init : vertical dynamics physics scheme' 118 WRITE(numout,*) '~~~~~~~~~~~' 119 IF( nzdf == 0 ) WRITE(numout,*) ' ===>> Explicit time-splitting scheme' 120 IF( nzdf == 1 ) WRITE(numout,*) ' ===>> Implicit (euler backward) scheme' 121 ENDIF 122 ! 123 END SUBROUTINE dyn_zdf_init 118 ! 119 IF( nn_timing == 1 ) CALL timing_start('dyn_zdf_imp') 120 ! 121 IF( kt == nit000 ) THEN 122 IF(lwp) WRITE(numout,*) 123 IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator' 124 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 125 ! 126 If( ln_linssh ) THEN ; r_vvl = 0._wp ! non-linear free surface indicator 127 ELSE ; r_vvl = 1._wp 128 ENDIF 129 ENDIF 130 ! 131 ! !== Time step dynamics ==! 132 ! 133 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity 134 DO jk = 1, jpkm1 135 ua(:,:,jk) = ( ub(:,:,jk) + p2dt * ua(:,:,jk) ) * umask(:,:,jk) 136 va(:,:,jk) = ( vb(:,:,jk) + p2dt * va(:,:,jk) ) * vmask(:,:,jk) 137 END DO 138 ELSE ! applied on thickness weighted velocity 139 DO jk = 1, jpkm1 140 ua(:,:,jk) = ( e3u_b(:,:,jk) * ub(:,:,jk) & 141 & + p2dt * e3u_n(:,:,jk) * ua(:,:,jk) ) / e3u_a(:,:,jk) * umask(:,:,jk) 142 va(:,:,jk) = ( e3v_b(:,:,jk) * vb(:,:,jk) & 143 & + p2dt * e3v_n(:,:,jk) * va(:,:,jk) ) / e3v_a(:,:,jk) * vmask(:,:,jk) 144 END DO 145 ENDIF 146 ! 147 ! !== Apply semi-implicit bottom friction ==! 148 ! 149 ! Only needed for semi-implicit bottom friction setup. The explicit 150 ! bottom friction has been included in "u(v)a" which act as the R.H.S 151 ! column vector of the tri-diagonal matrix equation 152 ! 153 IF( ln_drgimp ) THEN 154 DO jj = 2, jpjm1 155 DO ji = 2, jpim1 156 ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points 157 ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 158 avmu(ji,jj,ikbu+1) = -0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * e3uw_n(ji,jj,ikbu+1) 159 avmv(ji,jj,ikbv+1) = -0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * e3vw_n(ji,jj,ikbv+1) 160 END DO 161 END DO 162 IF ( ln_isfcav ) THEN 163 DO jj = 2, jpjm1 164 DO ji = 2, jpim1 165 ikbu = miku(ji,jj) ! ocean top level at u- and v-points 166 ikbv = mikv(ji,jj) ! (first wet ocean u- and v-points) 167 ! top Cd is masked (=0 outside cavities) no need of test on mik>=2 168 avmu(ji,jj,ikbu) = -0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * e3uw_n(ji,jj,ikbu) 169 avmv(ji,jj,ikbv) = -0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * e3vw_n(ji,jj,ikbv) 170 END DO 171 END DO 172 END IF 173 ENDIF 174 ! 175 ! With split-explicit free surface, barotropic stress is treated explicitly 176 ! Update velocities at the bottom. 177 ! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does 178 ! not lead to the effective stress seen over the whole barotropic loop. 179 ! G. Madec : in linear free surface, e3u_a = e3u_n = e3u_0, so systematic use of e3u_a 180 IF( ln_drgimp .AND. ln_dynspg_ts ) THEN 181 DO jk = 1, jpkm1 ! remove barotropic velocities 182 ua(:,:,jk) = ( ua(:,:,jk) - ua_b(:,:) ) * umask(:,:,jk) 183 va(:,:,jk) = ( va(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk) 184 END DO 185 DO jj = 2, jpjm1 ! Add bottom/top stress due to barotropic component only 186 DO ji = fs_2, fs_jpim1 ! vector opt. 187 ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points 188 ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 189 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,ikbu) + r_vvl * e3u_a(ji,jj,ikbu) 190 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikbv) + r_vvl * e3v_a(ji,jj,ikbv) 191 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua 192 va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va 193 END DO 194 END DO 195 IF( ln_isfcav ) THEN ! Ocean cavities (ISF) 196 DO jj = 2, jpjm1 197 DO ji = fs_2, fs_jpim1 ! vector opt. 198 ikbu = miku(ji,jj) ! top ocean level at u- and v-points 199 ikbv = mikv(ji,jj) ! (first wet ocean u- and v-points) 200 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,ikbu) + r_vvl * e3u_a(ji,jj,ikbu) 201 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikbv) + r_vvl * e3v_a(ji,jj,ikbv) 202 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + p2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua 203 va(ji,jj,ikbv) = va(ji,jj,ikbv) + p2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va 204 END DO 205 END DO 206 END IF 207 ENDIF 208 ! 209 ! !== Vertical diffusion on u ==! 210 ! 211 ! Matrix and second member construction 212 ! bottom boundary condition: both zwi and zws must be masked as avmu can take 213 ! non zero value at the ocean bottom depending on the bottom friction used. 214 ! 215 DO jk = 1, jpkm1 ! Matrix 216 DO jj = 2, jpjm1 217 DO ji = fs_2, fs_jpim1 ! vector opt. 218 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at T-point 219 zzwi = - p2dt * avmu(ji,jj,jk ) / ( ze3ua * e3uw_n(ji,jj,jk ) ) 220 zzws = - p2dt * avmu(ji,jj,jk+1) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) 221 zwi(ji,jj,jk) = zzwi * wumask(ji,jj,jk ) 222 zws(ji,jj,jk) = zzws * wumask(ji,jj,jk+1) 223 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 224 END DO 225 END DO 226 END DO 227 DO jj = 2, jpjm1 ! Surface boundary conditions 228 DO ji = fs_2, fs_jpim1 ! vector opt. 229 zwi(ji,jj,1) = 0._wp 230 zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 231 END DO 232 END DO 233 234 ! Matrix inversion starting from the first level 235 !----------------------------------------------------------------------- 236 ! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk ) 237 ! 238 ! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 ) 239 ! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 ) 240 ! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 ) 241 ! ( ... )( ... ) ( ... ) 242 ! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk ) 243 ! 244 ! m is decomposed in the product of an upper and a lower triangular matrix 245 ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 246 ! The solution (the after velocity) is in ua 247 !----------------------------------------------------------------------- 248 ! 249 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 250 DO jj = 2, jpjm1 251 DO ji = fs_2, fs_jpim1 ! vector opt. 252 zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 253 END DO 254 END DO 255 END DO 256 ! 257 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==! 258 DO ji = fs_2, fs_jpim1 ! vector opt. 259 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 260 ua(ji,jj,1) = ua(ji,jj,1) + p2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 261 & / ( ze3ua * rau0 ) * umask(ji,jj,1) 262 END DO 263 END DO 264 DO jk = 2, jpkm1 265 DO jj = 2, jpjm1 266 DO ji = fs_2, fs_jpim1 267 ua(ji,jj,jk) = ua(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 268 END DO 269 END DO 270 END DO 271 ! 272 DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk ==! 273 DO ji = fs_2, fs_jpim1 ! vector opt. 274 ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 275 END DO 276 END DO 277 DO jk = jpk-2, 1, -1 278 DO jj = 2, jpjm1 279 DO ji = fs_2, fs_jpim1 280 ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 281 END DO 282 END DO 283 END DO 284 ! 285 ! !== Vertical diffusion on v ==! 286 ! 287 ! Matrix and second member construction 288 ! bottom boundary condition: both zwi and zws must be masked as avmv can take 289 ! non zero value at the ocean bottom depending on the bottom friction used 290 ! 291 DO jk = 1, jpkm1 ! Matrix 292 DO jj = 2, jpjm1 293 DO ji = fs_2, fs_jpim1 ! vector opt. 294 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at T-point 295 zzwi = - p2dt * avmv (ji,jj,jk ) / ( ze3va * e3vw_n(ji,jj,jk ) ) 296 zzws = - p2dt * avmv (ji,jj,jk+1) / ( ze3va * e3vw_n(ji,jj,jk+1) ) 297 zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk ) 298 zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) 299 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 300 END DO 301 END DO 302 END DO 303 DO jj = 2, jpjm1 ! Surface boundary conditions 304 DO ji = fs_2, fs_jpim1 ! vector opt. 305 zwi(ji,jj,1) = 0._wp 306 zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 307 END DO 308 END DO 309 310 ! Matrix inversion 311 !----------------------------------------------------------------------- 312 ! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk ) 313 ! 314 ! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 ) 315 ! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 ) 316 ! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 ) 317 ! ( ... )( ... ) ( ... ) 318 ! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk ) 319 ! 320 ! m is decomposed in the product of an upper and lower triangular matrix 321 ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 322 ! The solution (after velocity) is in 2d array va 323 !----------------------------------------------------------------------- 324 ! 325 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 326 DO jj = 2, jpjm1 327 DO ji = fs_2, fs_jpim1 ! vector opt. 328 zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 329 END DO 330 END DO 331 END DO 332 ! 333 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==! 334 DO ji = fs_2, fs_jpim1 ! vector opt. 335 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 336 va(ji,jj,1) = va(ji,jj,1) + p2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 337 & / ( ze3va * rau0 ) * vmask(ji,jj,1) 338 END DO 339 END DO 340 DO jk = 2, jpkm1 341 DO jj = 2, jpjm1 342 DO ji = fs_2, fs_jpim1 ! vector opt. 343 va(ji,jj,jk) = va(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 344 END DO 345 END DO 346 END DO 347 ! 348 DO jj = 2, jpjm1 !== third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk ==! 349 DO ji = fs_2, fs_jpim1 ! vector opt. 350 va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 351 END DO 352 END DO 353 DO jk = jpk-2, 1, -1 354 DO jj = 2, jpjm1 355 DO ji = fs_2, fs_jpim1 356 va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 357 END DO 358 END DO 359 END DO 360 ! 361 IF( nn_timing == 1 ) CALL timing_stop('dyn_zdf_imp') 362 ! 363 END SUBROUTINE dyn_zdf_imp 124 364 125 365 !!==============================================================================
Note: See TracChangeset
for help on using the changeset viewer.
