- Timestamp:
- 2017-12-13T15:58:53+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf.F90
r7753 r9019 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 !! ----------------------------------------------------------------------9 10 !!---------------------------------------------------------------------- 11 !! dyn_zdf : Update the momentum trend with the vertical diffusion12 !! dyn_zdf _init : initializations of the vertical diffusion scheme8 !! 4.0 ! 2017-06 (G. Madec) remove the explicit time-stepping option + avm at t-point 9 !!---------------------------------------------------------------------- 10 11 !!---------------------------------------------------------------------- 12 !! dyn_zdf : compute the after velocity through implicit calculation of vertical mixing 13 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 dynzdf_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 21 USE dynldf ,ONLY: nldf, np_lap_i ! dynamics: type of lateral mixing 22 USE dynldf_iso,ONLY: akzu, akzv ! dynamics: vertical component of rotated lateral mixing 19 23 USE ldfdyn ! lateral diffusion: eddy viscosity coef. 20 24 USE trd_oce ! trends: ocean variables … … 24 28 USE lib_mpp ! MPP library 25 29 USE prtctl ! Print control 26 USE wrk_nemo ! Memory Allocation27 30 USE timing ! Timing 28 31 … … 30 33 PRIVATE 31 34 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 35 PUBLIC dyn_zdf ! routine called by step.F90 36 37 REAL(wp) :: r_vvl ! non-linear free surface indicator: =0 if ln_linssh=T, =1 otherwise 36 38 37 39 !! * Substitutions 38 40 # include "vectopt_loop_substitute.h90" 39 41 !!---------------------------------------------------------------------- 40 !! NEMO/OPA 3.3 , NEMO Consortium (2010)42 !! NEMO/OPA 4.0 , NEMO Consortium (2017) 41 43 !! $Id$ 42 44 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 43 45 !!---------------------------------------------------------------------- 44 45 46 CONTAINS 46 47 … … 49 50 !! *** ROUTINE dyn_zdf *** 50 51 !! 51 !! ** Purpose : compute the vertical ocean dynamics physics. 52 !! ** Purpose : compute the trend due to the vert. momentum diffusion 53 !! together with the Leap-Frog time stepping using an 54 !! implicit scheme. 55 !! 56 !! ** Method : - Leap-Frog time stepping on all trends but the vertical mixing 57 !! ua = ub + 2*dt * ua vector form or linear free surf. 58 !! ua = ( e3u_b*ub + 2*dt * e3u_n*ua ) / e3u_a otherwise 59 !! - update the after velocity with the implicit vertical mixing. 60 !! This requires to solver the following system: 61 !! ua = ua + 1/e3u_a dk+1[ mi(avm) / e3uw_a dk[ua] ] 62 !! with the following surface/top/bottom boundary condition: 63 !! surface: wind stress input (averaged over kt-1/2 & kt+1/2) 64 !! top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfdrg.F90) 65 !! 66 !! ** Action : (ua,va) after velocity 52 67 !!--------------------------------------------------------------------- 53 !! 54 INTEGER, INTENT( in ) :: kt ! ocean time-step index 55 ! 56 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 68 INTEGER, INTENT(in) :: kt ! ocean time-step index 69 ! 70 INTEGER :: ji, jj, jk ! dummy loop indices 71 INTEGER :: iku, ikv ! local integers 72 REAL(wp) :: zzwi, ze3ua, zdt ! local scalars 73 REAL(wp) :: zzws, ze3va ! - - 74 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi, zwd, zws ! 3D workspace 75 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdu, ztrdv ! - - 57 76 !!--------------------------------------------------------------------- 58 77 ! 59 IF( nn_timing == 1 ) CALL timing_start('dyn_zdf') 60 ! 61 ! ! set time step 78 IF( ln_timing ) CALL timing_start('dyn_zdf') 79 ! 80 IF( kt == nit000 ) THEN !* initialization 81 IF(lwp) WRITE(numout,*) 82 IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator' 83 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 84 ! 85 If( ln_linssh ) THEN ; r_vvl = 0._wp ! non-linear free surface indicator 86 ELSE ; r_vvl = 1._wp 87 ENDIF 88 ENDIF 89 ! !* set time step 62 90 IF( neuler == 0 .AND. kt == nit000 ) THEN ; r2dt = rdt ! = rdt (restart with Euler time stepping) 63 91 ELSEIF( kt <= nit000 + 1 ) THEN ; r2dt = 2. * rdt ! = 2 rdt (leapfrog) 64 92 ENDIF 65 93 66 IF( l_trddyn ) THEN !temporary save of ta and sa trends67 CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv)94 IF( l_trddyn ) THEN !* temporary save of ta and sa trends 95 ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) ) 68 96 ztrdu(:,:,:) = ua(:,:,:) 69 97 ztrdv(:,:,:) = va(:,:,:) 70 98 ENDIF 71 72 SELECT CASE ( nzdf ) ! compute lateral mixing trend and add it to the general trend 73 ! 74 CASE ( 0 ) ; CALL dyn_zdf_exp( kt, r2dt ) ! explicit scheme 75 CASE ( 1 ) ; CALL dyn_zdf_imp( kt, r2dt ) ! implicit scheme 76 ! 77 END SELECT 78 99 ! 100 ! !== RHS: Leap-Frog time stepping on all trends but the vertical mixing ==! (put in ua,va) 101 ! 102 ! ! time stepping except vertical diffusion 103 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity 104 DO jk = 1, jpkm1 105 ua(:,:,jk) = ( ub(:,:,jk) + r2dt * ua(:,:,jk) ) * umask(:,:,jk) 106 va(:,:,jk) = ( vb(:,:,jk) + r2dt * va(:,:,jk) ) * vmask(:,:,jk) 107 END DO 108 ELSE ! applied on thickness weighted velocity 109 DO jk = 1, jpkm1 110 ua(:,:,jk) = ( e3u_b(:,:,jk) * ub(:,:,jk) & 111 & + r2dt * e3u_n(:,:,jk) * ua(:,:,jk) ) / e3u_a(:,:,jk) * umask(:,:,jk) 112 va(:,:,jk) = ( e3v_b(:,:,jk) * vb(:,:,jk) & 113 & + r2dt * e3v_n(:,:,jk) * va(:,:,jk) ) / e3v_a(:,:,jk) * vmask(:,:,jk) 114 END DO 115 ENDIF 116 ! ! add top/bottom friction 117 ! With split-explicit free surface, barotropic stress is treated explicitly Update velocities at the bottom. 118 ! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does 119 ! not lead to the effective stress seen over the whole barotropic loop. 120 ! G. Madec : in linear free surface, e3u_a = e3u_n = e3u_0, so systematic use of e3u_a 121 IF( ln_drgimp .AND. ln_dynspg_ts ) THEN 122 DO jk = 1, jpkm1 ! remove barotropic velocities 123 ua(:,:,jk) = ( ua(:,:,jk) - ua_b(:,:) ) * umask(:,:,jk) 124 va(:,:,jk) = ( va(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk) 125 END DO 126 DO jj = 2, jpjm1 ! Add bottom/top stress due to barotropic component only 127 DO ji = fs_2, fs_jpim1 ! vector opt. 128 iku = mbku(ji,jj) ! ocean bottom level at u- and v-points 129 ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 130 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) 131 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) 132 ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua 133 va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va 134 END DO 135 END DO 136 IF( ln_isfcav ) THEN ! Ocean cavities (ISF) 137 DO jj = 2, jpjm1 138 DO ji = fs_2, fs_jpim1 ! vector opt. 139 iku = miku(ji,jj) ! top ocean level at u- and v-points 140 ikv = mikv(ji,jj) ! (first wet ocean u- and v-points) 141 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) 142 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) 143 ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua 144 va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va 145 END DO 146 END DO 147 END IF 148 ENDIF 149 ! 150 ! !== Vertical diffusion on u ==! 151 ! 152 ! !* Matrix construction 153 zdt = r2dt * 0.5 154 IF( nldf == np_lap_i ) THEN ! rotated lateral mixing: add its vertical mixing (akzu) 155 DO jk = 1, jpkm1 156 DO jj = 2, jpjm1 157 DO ji = fs_2, fs_jpim1 ! vector opt. 158 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at T-point 159 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) & 160 & / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk ) 161 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) & 162 & / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 163 zwi(ji,jj,jk) = zzwi 164 zws(ji,jj,jk) = zzws 165 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 166 END DO 167 END DO 168 END DO 169 ELSE ! standard case 170 DO jk = 1, jpkm1 171 DO jj = 2, jpjm1 172 DO ji = fs_2, fs_jpim1 ! vector opt. 173 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at T-point 174 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk ) 175 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 176 zwi(ji,jj,jk) = zzwi 177 zws(ji,jj,jk) = zzws 178 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 179 END DO 180 END DO 181 END DO 182 ENDIF 183 ! 184 DO jj = 2, jpjm1 !* Surface boundary conditions 185 DO ji = fs_2, fs_jpim1 ! vector opt. 186 zwi(ji,jj,1) = 0._wp 187 zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 188 END DO 189 END DO 190 ! 191 ! !== Apply semi-implicit bottom friction ==! 192 ! 193 ! Only needed for semi-implicit bottom friction setup. The explicit 194 ! bottom friction has been included in "u(v)a" which act as the R.H.S 195 ! column vector of the tri-diagonal matrix equation 196 ! 197 IF ( ln_drgimp ) THEN ! implicit bottom friction 198 DO jj = 2, jpjm1 199 DO ji = 2, jpim1 200 iku = mbku(ji,jj) ! ocean bottom level at u- and v-points 201 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) ! after scale factor at T-point 202 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua 203 END DO 204 END DO 205 IF ( ln_isfcav ) THEN ! top friction (always implicit) 206 DO jj = 2, jpjm1 207 DO ji = 2, jpim1 208 !!gm top Cd is masked (=0 outside cavities) no need of test on mik>=2 ==>> it has been suppressed 209 iku = miku(ji,jj) ! ocean top level at u- and v-points 210 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) ! after scale factor at T-point 211 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua 212 END DO 213 END DO 214 END IF 215 ENDIF 216 ! 217 ! Matrix inversion starting from the first level 218 !----------------------------------------------------------------------- 219 ! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk ) 220 ! 221 ! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 ) 222 ! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 ) 223 ! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 ) 224 ! ( ... )( ... ) ( ... ) 225 ! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk ) 226 ! 227 ! m is decomposed in the product of an upper and a lower triangular matrix 228 ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 229 ! The solution (the after velocity) is in ua 230 !----------------------------------------------------------------------- 231 ! 232 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 233 DO jj = 2, jpjm1 234 DO ji = fs_2, fs_jpim1 ! vector opt. 235 zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 236 END DO 237 END DO 238 END DO 239 ! 240 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==! 241 DO ji = fs_2, fs_jpim1 ! vector opt. 242 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1) 243 ua(ji,jj,1) = ua(ji,jj,1) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 244 & / ( ze3ua * rau0 ) * umask(ji,jj,1) 245 END DO 246 END DO 247 DO jk = 2, jpkm1 248 DO jj = 2, jpjm1 249 DO ji = fs_2, fs_jpim1 250 ua(ji,jj,jk) = ua(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 251 END DO 252 END DO 253 END DO 254 ! 255 DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk ==! 256 DO ji = fs_2, fs_jpim1 ! vector opt. 257 ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 258 END DO 259 END DO 260 DO jk = jpk-2, 1, -1 261 DO jj = 2, jpjm1 262 DO ji = fs_2, fs_jpim1 263 ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 264 END DO 265 END DO 266 END DO 267 ! 268 ! !== Vertical diffusion on v ==! 269 ! 270 ! !* Matrix construction 271 zdt = r2dt * 0.5 272 IF( nldf == np_lap_i ) THEN ! rotated lateral mixing: add its vertical mixing (akzu) 273 DO jk = 1, jpkm1 274 DO jj = 2, jpjm1 275 DO ji = fs_2, fs_jpim1 ! vector opt. 276 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at T-point 277 zzwi = - zdt * ( avm(ji,jj+1,jk )+ avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) & 278 & / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 279 zzws = - zdt * ( avm(ji,jj+1,jk+1)+ avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) & 280 & / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 281 zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk ) 282 zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) 283 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 284 END DO 285 END DO 286 END DO 287 ELSE ! standard case 288 DO jk = 1, jpkm1 289 DO jj = 2, jpjm1 290 DO ji = fs_2, fs_jpim1 ! vector opt. 291 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at T-point 292 zzwi = - zdt * ( avm(ji,jj+1,jk )+ avm(ji,jj,jk ) ) / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 293 zzws = - zdt * ( avm(ji,jj+1,jk+1)+ avm(ji,jj,jk+1) ) / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 294 zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk ) 295 zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) 296 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 297 END DO 298 END DO 299 END DO 300 ENDIF 301 ! 302 DO jj = 2, jpjm1 !* Surface boundary conditions 303 DO ji = fs_2, fs_jpim1 ! vector opt. 304 zwi(ji,jj,1) = 0._wp 305 zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 306 END DO 307 END DO 308 ! !== Apply semi-implicit top/bottom friction ==! 309 ! 310 ! Only needed for semi-implicit bottom friction setup. The explicit 311 ! bottom friction has been included in "u(v)a" which act as the R.H.S 312 ! column vector of the tri-diagonal matrix equation 313 ! 314 IF( ln_drgimp ) THEN 315 DO jj = 2, jpjm1 316 DO ji = 2, jpim1 317 ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 318 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) ! after scale factor at T-point 319 zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va 320 END DO 321 END DO 322 IF ( ln_isfcav ) THEN 323 DO jj = 2, jpjm1 324 DO ji = 2, jpim1 325 ikv = mikv(ji,jj) ! (first wet ocean u- and v-points) 326 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) ! after scale factor at T-point 327 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3va 328 END DO 329 END DO 330 ENDIF 331 ENDIF 332 333 ! Matrix inversion 334 !----------------------------------------------------------------------- 335 ! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk ) 336 ! 337 ! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 ) 338 ! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 ) 339 ! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 ) 340 ! ( ... )( ... ) ( ... ) 341 ! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk ) 342 ! 343 ! m is decomposed in the product of an upper and lower triangular matrix 344 ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 345 ! The solution (after velocity) is in 2d array va 346 !----------------------------------------------------------------------- 347 ! 348 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 349 DO jj = 2, jpjm1 350 DO ji = fs_2, fs_jpim1 ! vector opt. 351 zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 352 END DO 353 END DO 354 END DO 355 ! 356 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==! 357 DO ji = fs_2, fs_jpim1 ! vector opt. 358 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1) 359 va(ji,jj,1) = va(ji,jj,1) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 360 & / ( ze3va * rau0 ) * vmask(ji,jj,1) 361 END DO 362 END DO 363 DO jk = 2, jpkm1 364 DO jj = 2, jpjm1 365 DO ji = fs_2, fs_jpim1 ! vector opt. 366 va(ji,jj,jk) = va(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 367 END DO 368 END DO 369 END DO 370 ! 371 DO jj = 2, jpjm1 !== third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk ==! 372 DO ji = fs_2, fs_jpim1 ! vector opt. 373 va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 374 END DO 375 END DO 376 DO jk = jpk-2, 1, -1 377 DO jj = 2, jpjm1 378 DO ji = fs_2, fs_jpim1 379 va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 380 END DO 381 END DO 382 END DO 383 ! 79 384 IF( l_trddyn ) THEN ! save the vertical diffusive trends for further diagnostics 80 385 ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / r2dt - ztrdu(:,:,:) 81 386 ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / r2dt - ztrdv(:,:,:) 82 387 CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt ) 83 CALL wrk_dealloc( jpi, jpj, jpk,ztrdu, ztrdv )388 DEALLOCATE( ztrdu, ztrdv ) 84 389 ENDIF 85 390 ! ! print mean trends (used for debugging) … … 87 392 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 88 393 ! 89 IF( nn_timing == 1) CALL timing_stop('dyn_zdf')394 IF( ln_timing ) CALL timing_stop('dyn_zdf') 90 395 ! 91 396 END SUBROUTINE dyn_zdf 92 93 94 SUBROUTINE dyn_zdf_init95 !!----------------------------------------------------------------------96 !! *** ROUTINE dyn_zdf_init ***97 !!98 !! ** Purpose : initializations of the vertical diffusion scheme99 !!100 !! ** Method : implicit (euler backward) scheme (default)101 !! explicit (time-splitting) scheme if ln_zdfexp=T102 !!----------------------------------------------------------------------103 USE zdftke104 USE zdfgls105 !!----------------------------------------------------------------------106 !107 ! Choice from ln_zdfexp read in namelist in zdfini108 IF( ln_zdfexp ) THEN ; nzdf = 0 ! use explicit scheme109 ELSE ; nzdf = 1 ! use implicit scheme110 ENDIF111 !112 ! Force implicit schemes113 IF( lk_zdftke .OR. lk_zdfgls ) nzdf = 1 ! TKE or GLS physics114 IF( ln_dynldf_iso ) nzdf = 1 ! iso-neutral lateral physics115 IF( ln_dynldf_hor .AND. ln_sco ) nzdf = 1 ! horizontal lateral physics in s-coordinate116 !117 IF(lwp) THEN ! Print the choice118 WRITE(numout,*)119 WRITE(numout,*) 'dyn_zdf_init : vertical dynamics physics scheme'120 WRITE(numout,*) '~~~~~~~~~~~'121 IF( nzdf == 0 ) WRITE(numout,*) ' ===>> Explicit time-splitting scheme'122 IF( nzdf == 1 ) WRITE(numout,*) ' ===>> Implicit (euler backward) scheme'123 ENDIF124 !125 END SUBROUTINE dyn_zdf_init126 397 127 398 !!==============================================================================
Note: See TracChangeset
for help on using the changeset viewer.
