- Timestamp:
- 2020-07-03T20:46:53+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/OCE/DYN/dynzdf.F90
r12939 r13248 38 38 !! * Substitutions 39 39 # include "do_loop_substitute.h90" 40 # include "domzgr_substitute.h90" 40 41 !!---------------------------------------------------------------------- 41 42 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 55 56 !! ** Method : - Leap-Frog time stepping on all trends but the vertical mixing 56 57 !! u(after) = u(before) + 2*dt * u(rhs) vector form or linear free surf. 57 !! u(after) = ( e3u_b*u(before) + 2*dt * e3u_n*u(rhs) ) / e3u (after)otherwise58 !! u(after) = ( e3u_b*u(before) + 2*dt * e3u_n*u(rhs) ) / e3u_after otherwise 58 59 !! - update the after velocity with the implicit vertical mixing. 59 60 !! This requires to solver the following system: 60 !! u(after) = u(after) + 1/e3u (after) dk+1[ mi(avm) / e3uw(after)dk[ua] ]61 !! u(after) = u(after) + 1/e3u_after dk+1[ mi(avm) / e3uw_after dk[ua] ] 61 62 !! with the following surface/top/bottom boundary condition: 62 63 !! surface: wind stress input (averaged over kt-1/2 & kt+1/2) … … 112 113 ELSE ! applied on thickness weighted velocity 113 114 DO_3D_00_00( 1, jpkm1 ) 114 puu(ji,jj,jk,Kaa) = ( e3u(ji,jj,jk,Kbb) * puu(ji,jj,jk,Kbb) & 115 & + rDt * e3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Krhs) ) / e3u(ji,jj,jk,Kaa) * umask(ji,jj,jk) 116 pvv(ji,jj,jk,Kaa) = ( e3v(ji,jj,jk,Kbb) * pvv(ji,jj,jk,Kbb) & 117 & + rDt * e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Krhs) ) / e3v(ji,jj,jk,Kaa) * vmask(ji,jj,jk) 115 puu(ji,jj,jk,Kaa) = ( e3u(ji,jj,jk,Kbb) * puu(ji,jj,jk,Kbb ) & 116 & + rDt * e3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Krhs) ) & 117 & / e3u(ji,jj,jk,Kaa) * umask(ji,jj,jk) 118 pvv(ji,jj,jk,Kaa) = ( e3v(ji,jj,jk,Kbb) * pvv(ji,jj,jk,Kbb ) & 119 & + rDt * e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Krhs) ) & 120 & / e3v(ji,jj,jk,Kaa) * vmask(ji,jj,jk) 118 121 END_3D 119 122 ENDIF … … 131 134 iku = mbku(ji,jj) ! ocean bottom level at u- and v-points 132 135 ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 133 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa) 134 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa) 136 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) & 137 & + r_vvl * e3u(ji,jj,iku,Kaa) 138 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) & 139 & + r_vvl * e3v(ji,jj,ikv,Kaa) 135 140 puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + rDt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua 136 141 pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + rDt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va … … 140 145 iku = miku(ji,jj) ! top ocean level at u- and v-points 141 146 ikv = mikv(ji,jj) ! (first wet ocean u- and v-points) 142 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa) 143 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa) 147 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) & 148 & + r_vvl * e3u(ji,jj,iku,Kaa) 149 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) & 150 & + r_vvl * e3v(ji,jj,ikv,Kaa) 144 151 puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + rDt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua 145 152 pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + rDt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va … … 156 163 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu) 157 164 DO_3D_00_00( 1, jpkm1 ) 158 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point 165 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) & 166 & + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point 159 167 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) & 160 168 & / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) … … 169 177 CASE DEFAULT ! iso-level lateral mixing 170 178 DO_3D_00_00( 1, jpkm1 ) 171 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point 172 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) 173 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 179 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) & ! after scale factor at U-point 180 & + r_vvl * e3u(ji,jj,jk,Kaa) 181 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) & 182 & / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) 183 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) & 184 & / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 174 185 zWui = ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) / ze3ua 175 186 zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua … … 181 192 DO_2D_00_00 182 193 zwi(ji,jj,1) = 0._wp 183 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) + r_vvl * e3u(ji,jj,1,Kaa) 184 zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji ,jj,2) ) / ( ze3ua * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2) 194 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) & 195 & + r_vvl * e3u(ji,jj,1,Kaa) 196 zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji ,jj,2) ) & 197 & / ( ze3ua * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2) 185 198 zWus = ( wi(ji ,jj,2) + wi(ji+1,jj,2) ) / ze3ua 186 199 zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp ) … … 191 204 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu) 192 205 DO_3D_00_00( 1, jpkm1 ) 193 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point 206 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) & 207 & + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point 194 208 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) & 195 209 & / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) … … 202 216 CASE DEFAULT ! iso-level lateral mixing 203 217 DO_3D_00_00( 1, jpkm1 ) 204 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point 205 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) 206 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 218 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) & 219 & + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point 220 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) & 221 & / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) 222 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) & 223 & / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 207 224 zwi(ji,jj,jk) = zzwi 208 225 zws(ji,jj,jk) = zzws … … 226 243 DO_2D_00_00 227 244 iku = mbku(ji,jj) ! ocean bottom level at u- and v-points 228 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa) ! after scale factor at T-point 245 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) & 246 & + r_vvl * e3u(ji,jj,iku,Kaa) ! after scale factor at T-point 229 247 zwd(ji,jj,iku) = zwd(ji,jj,iku) - rDt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua 230 248 END_2D … … 233 251 !!gm top Cd is masked (=0 outside cavities) no need of test on mik>=2 ==>> it has been suppressed 234 252 iku = miku(ji,jj) ! ocean top level at u- and v-points 235 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa) ! after scale factor at T-point 253 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) & 254 & + r_vvl * e3u(ji,jj,iku,Kaa) ! after scale factor at T-point 236 255 zwd(ji,jj,iku) = zwd(ji,jj,iku) - rDt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua 237 256 END_2D … … 259 278 ! 260 279 DO_2D_00_00 261 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) + r_vvl * e3u(ji,jj,1,Kaa) 280 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) & 281 & + r_vvl * e3u(ji,jj,1,Kaa) 262 282 puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + rDt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 263 283 & / ( ze3ua * rho0 ) * umask(ji,jj,1) … … 282 302 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzv) 283 303 DO_3D_00_00( 1, jpkm1 ) 284 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 304 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) & 305 & + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 285 306 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) & 286 307 & / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk ) … … 295 316 CASE DEFAULT ! iso-level lateral mixing 296 317 DO_3D_00_00( 1, jpkm1 ) 297 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 298 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk ) 299 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 318 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) & 319 & + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 320 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) & 321 & / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk ) 322 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) & 323 & / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 300 324 zWvi = ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) / ze3va 301 325 zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va … … 307 331 DO_2D_00_00 308 332 zwi(ji,jj,1) = 0._wp 309 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) + r_vvl * e3v(ji,jj,1,Kaa) 310 zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2) 333 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) & 334 & + r_vvl * e3v(ji,jj,1,Kaa) 335 zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) & 336 & / ( ze3va * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2) 311 337 zWvs = ( wi(ji,jj ,2) + wi(ji,jj+1,2) ) / ze3va 312 338 zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp ) … … 317 343 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu) 318 344 DO_3D_00_00( 1, jpkm1 ) 319 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 345 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) & 346 & + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 320 347 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) & 321 348 & / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk ) … … 328 355 CASE DEFAULT ! iso-level lateral mixing 329 356 DO_3D_00_00( 1, jpkm1 ) 330 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 331 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk ) 332 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 357 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) & 358 & + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 359 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) & 360 & / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk ) 361 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) & 362 & / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 333 363 zwi(ji,jj,jk) = zzwi 334 364 zws(ji,jj,jk) = zzws … … 351 381 DO_2D_00_00 352 382 ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 353 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa) ! after scale factor at T-point 383 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) & 384 & + r_vvl * e3v(ji,jj,ikv,Kaa) ! after scale factor at T-point 354 385 zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - rDt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va 355 386 END_2D … … 357 388 DO_2D_00_00 358 389 ikv = mikv(ji,jj) ! (first wet ocean u- and v-points) 359 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa) ! after scale factor at T-point 390 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) & 391 & + r_vvl * e3v(ji,jj,ikv,Kaa) ! after scale factor at T-point 360 392 zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - rDt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / ze3va 361 393 END_2D … … 383 415 ! 384 416 DO_2D_00_00 385 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) + r_vvl * e3v(ji,jj,1,Kaa) 417 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) & 418 & + r_vvl * e3v(ji,jj,1,Kaa) 386 419 pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + rDt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 387 420 & / ( ze3va * rho0 ) * vmask(ji,jj,1)
Note: See TracChangeset
for help on using the changeset viewer.
