Changeset 9976 for NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP
- Timestamp:
- 2018-07-19T17:58:12+02:00 (6 years ago)
- Location:
- NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP
- Files:
-
- 12 edited
-
cfgs/SHARED/field_def_nemo-oce.xml (modified) (1 diff)
-
cfgs/SHARED/namelist_ref (modified) (1 diff)
-
src/OCE/DOM/dom_oce.F90 (modified) (2 diffs)
-
src/OCE/DOM/domvvl.F90 (modified) (2 diffs)
-
src/OCE/DYN/dynzdf.F90 (modified) (3 diffs)
-
src/OCE/DYN/sshwzv.F90 (modified) (3 diffs)
-
src/OCE/TRA/trazdf.F90 (modified) (1 diff)
-
src/OCE/ZDF/zdf_oce.F90 (modified) (1 diff)
-
src/OCE/ZDF/zdfosm.F90 (modified) (1 diff)
-
src/OCE/ZDF/zdfphy.F90 (modified) (4 diffs)
-
src/OCE/oce.F90 (modified) (2 diffs)
-
src/OCE/step.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/cfgs/SHARED/field_def_nemo-oce.xml
r9957 r9976 78 78 <field id="Lf_NHpf" long_name="MLE: Lf = N H / f" unit="m" /> 79 79 80 <!-- next variables available with ln_zad_Aimp=.true. --> 81 <field id="Courant" long_name="Courant number" unit="#" grid_ref="grid_T_3D" /> 82 <field id="wimp" long_name="Implicit vertical velocity" unit="m/s" grid_ref="grid_T_3D" /> 83 <field id="wexp" long_name="Explicit vertical velocity" unit="m/s" grid_ref="grid_T_3D" /> 80 84 <!-- next variables available with key_diahth --> 81 85 <field id="mlddzt" long_name="Thermocline Depth (depth of max dT/dz)" standard_name="depth_at_maximum_upward_derivative_of_sea_water_potential_temperature" unit="m" /> -
NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/cfgs/SHARED/namelist_ref
r9957 r9976 972 972 &namzdf ! vertical physics manager (default: NO selection) 973 973 !----------------------------------------------------------------------- 974 ! ! adaptive-implicit vertical advection 975 ln_zad_Aimp = .false. ! Courant number dependent scheme (Shchepetkin 2015) 976 ! 974 977 ! ! type of vertical closure (required) 975 978 ln_zdfcst = .false. ! constant mixing -
NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/DOM/dom_oce.F90
r9957 r9976 135 135 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3f_0 , e3f_n !: f- vert. scale factor [m] 136 136 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3w_0 , e3w_b , e3w_n !: w- vert. scale factor [m] 137 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3uw_0 , e3uw_b , e3uw_n 138 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3vw_0 , e3vw_b , e3vw_n 137 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3uw_0 , e3uw_b , e3uw_n , e3uw_a !: uw-vert. scale factor [m] 138 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e3vw_0 , e3vw_b , e3vw_n , e3vw_a !: vw-vert. scale factor [m] 139 139 140 140 ! ! ref. ! before ! now ! … … 267 267 & e3uw_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) , & 268 268 & e3uw_b(jpi,jpj,jpk) , e3vw_b(jpi,jpj,jpk) , & 269 & e3uw_n(jpi,jpj,jpk) , e3vw_n(jpi,jpj,jpk) , STAT=ierr(5) ) 269 & e3uw_n(jpi,jpj,jpk) , e3vw_n(jpi,jpj,jpk) , & 270 & e3uw_a(jpi,jpj,jpk) , e3vw_a(jpi,jpj,jpk) , STAT=ierr(5) ) 270 271 ! 271 272 ALLOCATE( ht_0(jpi,jpj) , hu_0(jpi,jpj) , hv_0(jpi,jpj) , & -
NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/DOM/domvvl.F90
r9957 r9976 149 149 150 150 ! We need to define e3[tuv]_a for AGRIF initialisation (should not be a problem for the restartability...) 151 e3t_a(:,:,:) = e3t_n(:,:,:) 152 e3u_a(:,:,:) = e3u_n(:,:,:) 153 e3v_a(:,:,:) = e3v_n(:,:,:) 151 e3t_a (:,:,:) = e3t_n (:,:,:) 152 e3u_a (:,:,:) = e3u_n (:,:,:) 153 e3v_a (:,:,:) = e3v_n (:,:,:) 154 e3uw_a(:,:,:) = e3uw_n(:,:,:) 155 e3vw_a(:,:,:) = e3vw_n(:,:,:) 154 156 ! 155 157 ! !== depth of t and w-point ==! (set the isf depth as it is in the initial timestep) … … 534 536 END IF 535 537 536 ! *********************************** !537 ! After scale factors at u- v- points !538 ! *********************************** !538 ! ******************************************* ! 539 ! After scale factors at u- v- uw- vw- points ! 540 ! ******************************************* ! 539 541 540 542 CALL dom_vvl_interpol( e3t_a(:,:,:), e3u_a(:,:,:), 'U' ) 541 543 CALL dom_vvl_interpol( e3t_a(:,:,:), e3v_a(:,:,:), 'V' ) 544 CALL dom_vvl_interpol( e3u_a(:,:,:), e3uw_a(:,:,:), 'UW' ) 545 CALL dom_vvl_interpol( e3v_a(:,:,:), e3vw_a(:,:,:), 'VW' ) 542 546 543 547 ! *********************************** ! -
NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/DYN/dynzdf.F90
r9957 r9976 71 71 REAL(wp) :: zzwi, ze3ua, zdt ! local scalars 72 72 REAL(wp) :: zzws, ze3va ! - - 73 REAL(wp) :: z1_e3un, z1_e3vn ! - - 74 REAL(wp) :: zWu , zWv ! - - 75 REAL(wp) :: zWui, zWvi ! - - 76 REAL(wp) :: zWus, zWvs ! - - 73 77 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi, zwd, zws ! 3D workspace 74 78 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdu, ztrdv ! - - … … 155 159 ! !* Matrix construction 156 160 zdt = r2dt * 0.5 157 SELECT CASE( nldf_dyn ) 158 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu) 159 DO jk = 1, jpkm1 160 DO jj = 2, jpjm1 161 DO ji = fs_2, fs_jpim1 ! vector opt. 162 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at T-point 163 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) & 164 & / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk ) 165 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) & 166 & / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 167 zwi(ji,jj,jk) = zzwi 168 zws(ji,jj,jk) = zzws 169 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 170 END DO 171 END DO 172 END DO 173 CASE DEFAULT ! iso-level lateral mixing 174 DO jk = 1, jpkm1 175 DO jj = 2, jpjm1 176 DO ji = fs_2, fs_jpim1 ! vector opt. 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 = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk ) 179 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) 180 zwi(ji,jj,jk) = zzwi 181 zws(ji,jj,jk) = zzws 182 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 183 END DO 184 END DO 185 END DO 186 END SELECT 161 IF( ln_zad_Aimp ) THEN !! 162 IF( ln_dynadv_vec ) THEN !== Vector invariant advection ==! 163 SELECT CASE( nldf_dyn ) 164 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu) 165 DO jk = 1, jpkm1 166 DO jj = 2, jpjm1 167 DO ji = fs_2, fs_jpim1 ! vector opt. 168 z1_e3un = 1._wp / e3u_n(ji,jj,jk) 169 zzwi = ( ( avm (ji+1,jj,jk ) + avm (ji,jj,jk ) + akzu(ji,jj,jk ) ) & 170 & / e3uw_a(ji ,jj,jk ) ) * z1_e3un * wumask(ji,jj,jk ) 171 zzws = ( ( avm (ji+1,jj,jk+1) + avm (ji,jj,jk+1) + akzu(ji,jj,jk+1) ) & 172 & / e3uw_a(ji ,jj,jk+1) ) * z1_e3un * wumask(ji,jj,jk+1) 173 zWu = 0.25_wp * ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) & 174 & + wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) 175 zwi(ji,jj,jk) = - zdt * ( zzwi + MIN( zWu, 0._wp ) * z1_e3un ) 176 zws(ji,jj,jk) = - zdt * ( zzws - MAX( zWu, 0._wp ) * z1_e3un ) 177 zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws + ( - MAX( zWu, 0._wp ) + MIN( zWu, 0._wp ) ) * z1_e3un ) 178 END DO 179 END DO 180 END DO 181 CASE DEFAULT ! iso-level lateral mixing 182 DO jk = 1, jpkm1 183 DO jj = 2, jpjm1 184 DO ji = fs_2, fs_jpim1 ! vector opt. 185 z1_e3un = 1._wp / e3u_n(ji,jj,jk) 186 zzwi = ( ( avm (ji+1,jj,jk ) + avm(ji,jj,jk ) ) & 187 & / e3uw_a(ji ,jj,jk ) ) * z1_e3un * wumask(ji,jj,jk ) 188 zzws = ( ( avm (ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) & 189 & / e3uw_a(ji ,jj,jk+1) ) * z1_e3un * wumask(ji,jj,jk+1) 190 zWu = 0.25_wp * ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) & 191 & + wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) 192 zwi(ji,jj,jk) = - zdt * ( zzwi + MIN( zWu, 0._wp ) * z1_e3un ) 193 zws(ji,jj,jk) = - zdt * ( zzws - MAX( zWu, 0._wp ) * z1_e3un ) 194 zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws + ( - MAX( zWu, 0._wp ) + MIN( zWu, 0._wp ) ) * z1_e3un ) 195 END DO 196 END DO 197 END DO 198 END SELECT 199 ELSE !== Flux form advection ==! 200 SELECT CASE( nldf_dyn ) 201 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu) 202 DO jk = 1, jpkm1 203 DO jj = 2, jpjm1 204 DO ji = fs_2, fs_jpim1 ! vector opt. 205 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point 206 zzwi = ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) & 207 & / ( ze3ua * e3uw_a(ji,jj,jk ) ) * wumask(ji,jj,jk ) 208 zzws = ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) & 209 & / ( ze3ua * e3uw_a(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 210 zWui = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) 211 zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) 212 zwi(ji,jj,jk) = - zdt * ( zzwi + MIN( zWui, 0._wp ) ) 213 zws(ji,jj,jk) = - zdt * ( zzws - MAX( zWus, 0._wp ) ) 214 zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws - MAX( zWui, 0._wp ) + MIN( zWus, 0._wp ) ) 215 END DO 216 END DO 217 END DO 218 CASE DEFAULT ! iso-level lateral mixing 219 DO jk = 1, jpkm1 220 DO jj = 2, jpjm1 221 DO ji = fs_2, fs_jpim1 ! vector opt. 222 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point 223 zzwi = ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw_a(ji,jj,jk ) ) * wumask(ji,jj,jk ) 224 zzws = ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_a(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 225 zWui = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) 226 zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) 227 zwi(ji,jj,jk) = - zdt * ( zzwi + MIN( zWui, 0._wp ) ) 228 zws(ji,jj,jk) = - zdt * ( zzws - MAX( zWus, 0._wp ) ) 229 zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws - MAX( zWui, 0._wp ) + MIN( zWus, 0._wp ) ) 230 END DO 231 END DO 232 END DO 233 END SELECT 234 ENDIF 235 ELSE 236 SELECT CASE( nldf_dyn ) 237 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu) 238 DO jk = 1, jpkm1 239 DO jj = 2, jpjm1 240 DO ji = fs_2, fs_jpim1 ! vector opt. 241 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point 242 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) & 243 & / ( ze3ua * e3uw_a(ji,jj,jk ) ) * wumask(ji,jj,jk ) 244 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) & 245 & / ( ze3ua * e3uw_a(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 246 zwi(ji,jj,jk) = zzwi 247 zws(ji,jj,jk) = zzws 248 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 249 END DO 250 END DO 251 END DO 252 CASE DEFAULT ! iso-level lateral mixing 253 DO jk = 1, jpkm1 254 DO jj = 2, jpjm1 255 DO ji = fs_2, fs_jpim1 ! vector opt. 256 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point 257 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw_a(ji,jj,jk ) ) * wumask(ji,jj,jk ) 258 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw_a(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 259 zwi(ji,jj,jk) = zzwi 260 zws(ji,jj,jk) = zzws 261 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 262 END DO 263 END DO 264 END DO 265 END SELECT 266 ENDIF 187 267 ! 188 268 DO jj = 2, jpjm1 !* Surface boundary conditions … … 274 354 ! !* Matrix construction 275 355 zdt = r2dt * 0.5 276 SELECT CASE( nldf_dyn ) 277 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu) 278 DO jk = 1, jpkm1 279 DO jj = 2, jpjm1 280 DO ji = fs_2, fs_jpim1 ! vector opt. 281 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at T-point 282 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) & 283 & / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 284 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) & 285 & / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 286 zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk ) 287 zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) 288 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 289 END DO 290 END DO 291 END DO 292 CASE DEFAULT ! iso-level lateral mixing 293 DO jk = 1, jpkm1 294 DO jj = 2, jpjm1 295 DO ji = fs_2, fs_jpim1 ! vector opt. 296 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at T-point 297 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 298 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) 299 zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk ) 300 zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) 301 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 302 END DO 303 END DO 304 END DO 305 END SELECT 356 IF( ln_zad_Aimp ) THEN !! 357 IF( ln_dynadv_vec ) THEN !== Vector invariant advection ==! 358 SELECT CASE( nldf_dyn ) 359 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzv) 360 DO jk = 1, jpkm1 361 DO jj = 2, jpjm1 362 DO ji = fs_2, fs_jpim1 ! vector opt. 363 z1_e3vn = 1._wp / e3v_n(ji,jj,jk) 364 zzwi = ( ( avm (ji,jj+1,jk ) + avm (ji,jj,jk ) + akzv(ji,jj,jk ) ) & 365 & / e3vw_a(ji,jj ,jk ) ) * z1_e3vn * wvmask(ji,jj,jk ) 366 zzws = ( ( avm (ji,jj+1,jk+1) + avm (ji,jj,jk+1) + akzv(ji,jj,jk+1) ) & 367 & / e3vw_a(ji,jj ,jk+1) ) * z1_e3vn * wvmask(ji,jj,jk+1) 368 zWv = 0.25_wp * ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) & 369 & + wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) 370 zwi(ji,jj,jk) = - zdt * ( zzwi + MIN( zWv, 0._wp ) * z1_e3vn ) 371 zws(ji,jj,jk) = - zdt * ( zzws - MAX( zWv, 0._wp ) * z1_e3vn ) 372 zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws + ( - MAX( zWv, 0._wp ) + MIN( zWv, 0._wp ) ) * z1_e3vn ) 373 END DO 374 END DO 375 END DO 376 CASE DEFAULT ! iso-level lateral mixing 377 DO jk = 1, jpkm1 378 DO jj = 2, jpjm1 379 DO ji = fs_2, fs_jpim1 ! vector opt. 380 z1_e3vn = 1._wp / e3v_n(ji,jj,jk) 381 zzwi = ( ( avm (ji,jj+1,jk ) + avm(ji,jj,jk ) ) & 382 & / e3vw_a(ji,jj ,jk ) ) * z1_e3vn * wvmask(ji,jj,jk ) 383 zzws = ( ( avm (ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) & 384 & / e3vw_a(ji ,jj,jk+1) ) * z1_e3vn * wvmask(ji,jj,jk+1) 385 zWv = 0.25_wp * ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) & 386 & + wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) 387 zwi(ji,jj,jk) = - zdt * ( zzwi + MIN( zWv, 0._wp ) * z1_e3vn ) 388 zws(ji,jj,jk) = - zdt * ( zzws - MAX( zWv, 0._wp ) * z1_e3vn ) 389 zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws + ( - MAX( zWv, 0._wp ) + MIN( zWv, 0._wp ) ) * z1_e3vn ) 390 END DO 391 END DO 392 END DO 393 END SELECT 394 ELSE !== Flux form advection ==! 395 SELECT CASE( nldf_dyn ) 396 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzv) 397 DO jk = 1, jpkm1 398 DO jj = 2, jpjm1 399 DO ji = fs_2, fs_jpim1 ! vector opt. 400 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at U-point 401 zzwi = ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) & 402 & / ( ze3va * e3vw_a(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 403 zzws = ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) & 404 & / ( ze3va * e3vw_a(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 405 zWvi = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) 406 zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) 407 zwi(ji,jj,jk) = - zdt * ( zzwi + MIN( zWvi, 0._wp ) ) 408 zws(ji,jj,jk) = - zdt * ( zzws - MAX( zWvs, 0._wp ) ) 409 zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) ) 410 END DO 411 END DO 412 END DO 413 CASE DEFAULT ! iso-level lateral mixing 414 DO jk = 1, jpkm1 415 DO jj = 2, jpjm1 416 DO ji = fs_2, fs_jpim1 ! vector opt. 417 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at U-point 418 zzwi = ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw_a(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 419 zzws = ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw_a(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 420 zWvi = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) 421 zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) 422 zwi(ji,jj,jk) = - zdt * ( zzwi + MIN( zWvi, 0._wp ) ) 423 zws(ji,jj,jk) = - zdt * ( zzws - MAX( zWvs, 0._wp ) ) 424 zwd(ji,jj,jk) = 1._wp + zdt * ( zzwi + zzws - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) ) 425 END DO 426 END DO 427 END DO 428 END SELECT 429 ENDIF 430 ELSE 431 SELECT CASE( nldf_dyn ) 432 CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu) 433 DO jk = 1, jpkm1 434 DO jj = 2, jpjm1 435 DO ji = fs_2, fs_jpim1 ! vector opt. 436 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at T-point 437 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) & 438 & / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 439 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) & 440 & / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 441 zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk ) 442 zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) 443 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 444 END DO 445 END DO 446 END DO 447 CASE DEFAULT ! iso-level lateral mixing 448 DO jk = 1, jpkm1 449 DO jj = 2, jpjm1 450 DO ji = fs_2, fs_jpim1 ! vector opt. 451 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at T-point 452 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 453 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) 454 zwi(ji,jj,jk) = zzwi * wvmask(ji,jj,jk ) 455 zws(ji,jj,jk) = zzws * wvmask(ji,jj,jk+1) 456 zwd(ji,jj,jk) = 1._wp - zzwi - zzws 457 END DO 458 END DO 459 END DO 460 END SELECT 461 ENDIF 306 462 ! 307 463 DO jj = 2, jpjm1 !* Surface boundary conditions -
NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/DYN/sshwzv.F90
r9957 r9976 28 28 #endif 29 29 ! 30 USE iom 30 31 USE in_out_manager ! I/O manager 31 32 USE restart ! only for lrst_oce … … 41 42 PUBLIC ssh_nxt ! called by step.F90 42 43 PUBLIC wzv ! called by step.F90 44 PUBLIC wAimp ! called by step.F90 43 45 PUBLIC ssh_swp ! called by step.F90 44 46 … … 265 267 END SUBROUTINE ssh_swp 266 268 269 SUBROUTINE wAimp( kt ) 270 !!---------------------------------------------------------------------- 271 !! *** ROUTINE wAimp *** 272 !! 273 !! ** Purpose : compute the Courant number and partition vertical velocity 274 !! if a proportion needs to be treated implicitly 275 !! 276 !! ** Method : - 277 !! 278 !! ** action : wn : now vertical velocity (to be handled explicitly) 279 !! : wi : now vertical velocity (for implicit treatment) 280 !! 281 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 282 !!---------------------------------------------------------------------- 283 INTEGER, INTENT(in) :: kt ! time step 284 ! 285 INTEGER :: ji, jj, jk ! dummy loop indices 286 REAL(wp) :: zCu, zcff, z1_e3w ! local scalars 287 REAL(wp) :: zcff_min, zcff_max ! local scalars 288 REAL(wp) , PARAMETER :: Cu_min = 0.6_wp ! local parameters 289 REAL(wp) , PARAMETER :: Cu_max = 0.95_wp ! local parameters 290 REAL(wp) , PARAMETER :: Cu_cut = 2._wp*Cu_max - Cu_min ! local parameters 291 REAL(wp) , PARAMETER :: Fcu = 4._wp*Cu_max*(Cu_max-Cu_min) ! local parameters 292 !!---------------------------------------------------------------------- 293 ! 294 IF( ln_timing ) CALL timing_start('wAimp') 295 ! 296 IF( kt == nit000 ) THEN 297 IF(lwp) WRITE(numout,*) 298 IF(lwp) WRITE(numout,*) 'wAimp : Courant number-based partitioning of now vertical velocity ' 299 IF(lwp) WRITE(numout,*) '~~~~~ ' 300 ! 301 wi (:,:,jpk) = 0._wp ! bottom value : wi=0 (set once for all) 302 Cu_adv(:,:,jpk) = 0._wp ! bottom value : Cu_adv=0 (set once for all) 303 ENDIF 304 ! 305 DO jk = 1, jpkm1 ! calculate Courant numbers 306 DO jj = 2, jpjm1 307 DO ji = 2, fs_jpim1 ! vector opt. 308 z1_e3w = 1._wp / e3w_n(ji,jj,jk) 309 Cu_adv(ji,jj,jk) = r2dt * ( ( MAX( wn(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) ) & 310 & + ( MAX( e2u(ji ,jj)*e3uw_n(ji ,jj,jk)*un(ji ,jj,jk), 0._wp ) - & 311 & MIN( e2u(ji-1,jj)*e3uw_n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) ) & 312 & / e2t(ji,jj) & 313 & + ( MAX( e1v(ji,jj )*e3vw_n(ji,jj ,jk)*vn(ji,jj ,jk), 0._wp ) - & 314 & MIN( e1v(ji,jj-1)*e3vw_n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) ) & 315 & / e1t(ji,jj) & 316 & ) * z1_e3w 317 END DO 318 END DO 319 END DO 320 ! 321 CALL iom_put("Courant",Cu_adv) 322 ! 323 wi(:,:,:) = 0._wp 324 zcff_min=999. ; zcff_max = -999. 325 IF( MAXVAL( Cu_adv(:,:,:) ) > Cu_min ) THEN ! Quick check if any breaches anywhere 326 DO jk = 1, jpkm1 ! or scan Courant criterion and partition 327 DO jj = 2, jpjm1 ! w where necessary 328 DO ji = 2, fs_jpim1 ! vector opt. 329 ! 330 zCu = MAX( Cu_adv(ji,jj,jk) , Cu_adv(ji,jj,jk+1) ) 331 ! 332 IF( zCu < Cu_min ) THEN !<-- Fully explicit 333 zcff = 0._wp 334 ELSEIF( zCu < Cu_cut ) THEN !<-- Mixed explicit 335 zcff = ( zCu - Cu_min )**2 336 zcff = zcff / ( Fcu + zcff ) 337 ELSE !<-- Fully implicit 338 zcff = ( zCu - Cu_max )/ zCu 339 ENDIF 340 ! 341 !zcff = 0._wp 342 wi(ji,jj,jk) = zcff * wn(ji,jj,jk) 343 wn(ji,jj,jk) = ( 1._wp - zcff ) * wn(ji,jj,jk) 344 zcff_min = MIN( zcff, zcff_min ) 345 zcff_max = MAX( zcff, zcff_max ) 346 ! 347 END DO 348 END DO 349 END DO 350 ELSE 351 ! Fully explicit everywhere 352 write(*,*) 'NoWi ',kt 353 ENDIF 354 write(*,*) 'ZCFF ',kt,zcff_min,zcff_max 355 CALL iom_put("wimp",wi) 356 CALL iom_put("wexp",wn) 357 ! 358 IF( ln_timing ) CALL timing_stop('wAimp') 359 ! 360 END SUBROUTINE wAimp 267 361 !!====================================================================== 268 362 END MODULE sshwzv -
NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/TRA/trazdf.F90
r9957 r9976 177 177 ! 178 178 ! Diagonal, lower (i), upper (s) (including the bottom boundary condition since avt is masked) 179 DO jk = 1, jpkm1 180 DO jj = 2, jpjm1 181 DO ji = fs_2, fs_jpim1 ! vector opt. 182 !!gm BUG I think, use e3w_a instead of e3w_n, not sure of that 183 zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk ) / e3w_n(ji,jj,jk ) 184 zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1) 185 zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zwi(ji,jj,jk) - zws(ji,jj,jk) 186 END DO 187 END DO 188 END DO 179 IF( ln_zad_Aimp ) THEN ! Adaptive implicit vertical advection 180 DO jk = 1, jpkm1 181 DO jj = 2, jpjm1 182 DO ji = fs_2, fs_jpim1 ! vector opt. 183 zwi(ji,jj,jk) = - p2dt * ( zwt(ji,jj,jk ) / e3w_n(ji,jj,jk ) + MIN( wi(ji,jj,jk ) , 0._wp ) ) 184 zws(ji,jj,jk) = - p2dt * ( zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1) - MAX( wi(ji,jj,jk+1) , 0._wp ) ) 185 zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zwi(ji,jj,jk) - zws(ji,jj,jk) & 186 & - p2dt * ( MAX( wi(ji,jj,jk ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) 187 END DO 188 END DO 189 END DO 190 ELSE 191 DO jk = 1, jpkm1 192 DO jj = 2, jpjm1 193 DO ji = fs_2, fs_jpim1 ! vector opt. 194 zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk ) / e3w_n(ji,jj,jk) 195 zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w_n(ji,jj,jk+1) 196 zwd(ji,jj,jk) = e3t_a(ji,jj,jk) - zwi(ji,jj,jk) - zws(ji,jj,jk) 197 END DO 198 END DO 199 END DO 200 ENDIF 189 201 ! 190 202 !! Matrix inversion from the first level -
NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/ZDF/zdf_oce.F90
r9957 r9976 18 18 19 19 ! !!* namelist namzdf: vertical physics * 20 ! ! Adaptive-implicit vertical advection flag 21 LOGICAL , PUBLIC :: ln_zad_Aimp !: adaptive (Courant number-based) implicit vertical advection 20 22 ! ! vertical closure scheme flags 21 23 LOGICAL , PUBLIC :: ln_zdfcst !: constant coefficients -
NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/ZDF/zdfosm.F90
r9957 r9976 32 32 33 33 !!---------------------------------------------------------------------- 34 !! ' key_zdfosm' OSMOSIS scheme34 !! 'ln_zdfosm' OSMOSIS scheme 35 35 !!---------------------------------------------------------------------- 36 36 !! zdf_osm : update momentum and tracer Kz from osm scheme -
NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/ZDF/zdfphy.F90
r9957 r9976 80 80 & ln_zdfswm, & ! surface wave-induced mixing 81 81 & ln_zdfiwm, & ! internal - - - 82 & ln_zad_Aimp, & ! apdative-implicit vertical advection 82 83 & rn_avm0, rn_avt0, nn_avb, nn_havtb ! coefficients 83 84 !!---------------------------------------------------------------------- … … 101 102 IF(lwp) THEN ! Parameter print 102 103 WRITE(numout,*) ' Namelist namzdf : set vertical mixing mixing parameters' 104 WRITE(numout,*) ' adaptive-implicit vertical advection' 105 WRITE(numout,*) ' Courant number targeted application ln_zad_Aimp = ', ln_zad_Aimp 103 106 WRITE(numout,*) ' vertical closure scheme' 104 107 WRITE(numout,*) ' constant vertical mixing coefficient ln_zdfcst = ', ln_zdfcst … … 127 130 ENDIF 128 131 132 IF( ln_zad_Aimp ) THEN 133 IF( zdf_phy_alloc() /= 0 ) & 134 & CALL ctl_stop( 'STOP', 'zdf_phy_init : unable to allocate adaptive-implicit z-advection arrays' ) 135 wi(:,:,:) = 0._wp 136 ENDIF 129 137 ! !== Background eddy viscosity and diffusivity ==! 130 138 IF( nn_avb == 0 ) THEN ! Define avmb, avtb from namelist parameter … … 316 324 ! 317 325 END SUBROUTINE zdf_phy 326 INTEGER FUNCTION zdf_phy_alloc() 327 !!---------------------------------------------------------------------- 328 !! *** FUNCTION zdf_phy_alloc *** 329 !!---------------------------------------------------------------------- 330 ! Allocate wi array (declared in oce.F90) for use with the adaptive-implicit vertical velocity option 331 ALLOCATE( wi(jpi,jpj,jpk), Cu_adv(jpi,jpj,jpk), STAT= zdf_phy_alloc ) 332 IF( zdf_phy_alloc /= 0 ) CALL ctl_warn('zdf_phy_alloc: failed to allocate ln_zad_Aimp=T required arrays') 333 IF( lk_mpp ) CALL mpp_sum ( zdf_phy_alloc ) 334 END FUNCTION zdf_phy_alloc 318 335 319 336 !!====================================================================== -
NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/oce.F90
r9957 r9976 22 22 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: vb , vn , va !: j-horizontal velocity [m/s] 23 23 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: wn !: vertical velocity [m/s] 24 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: wi !: vertical vel. (adaptive-implicit) [m/s] 24 25 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdivn !: horizontal divergence [s-1] 25 26 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: tsb , tsn , tsa !: 4D T-S fields [Celsius,psu] … … 29 30 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: rhd !: in situ density anomalie rhd=(rho-rau0)/rau0 [no units] 30 31 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: rhop !: potential volumic mass [kg/m3] 32 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: Cu_adv !: vertical Courant number (adaptive-implicit) 31 33 32 34 !! free surface ! before ! now ! after ! -
NEMO/branches/2018/dev_r9956_ENHANCE05_ZAD_AIMP/src/OCE/step.F90
r9957 r9976 160 160 161 161 CALL ssh_nxt ( kstp ) ! after ssh (includes call to div_hor) 162 IF( .NOT.ln_linssh )CALL dom_vvl_sf_nxt( kstp ) ! after vertical scale factors162 IF( .NOT.ln_linssh ) CALL dom_vvl_sf_nxt( kstp ) ! after vertical scale factors 163 163 CALL wzv ( kstp ) ! now cross-level velocity 164 IF( ln_zad_Aimp ) CALL wAimp ( kstp ) ! Adaptive-implicit vertical advection partitioning 164 165 CALL eos ( tsn, rhd, rhop, gdept_n(:,:,:) ) ! now in situ density for hpg computation 165 166 … … 198 199 IF(.NOT.ln_linssh) CALL dom_vvl_sf_nxt( kstp, kcall=2 ) ! after vertical scale factors (update depth average component) 199 200 CALL wzv ( kstp ) ! now cross-level velocity 201 IF( ln_zad_Aimp ) CALL wAimp ( kstp ) ! Adaptive-implicit vertical advection partitioning 200 202 ENDIF 201 203
Note: See TracChangeset
for help on using the changeset viewer.
