- Timestamp:
- 2018-07-23T11:33:03+02:00 (6 years ago)
- Location:
- branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/ZDF
- Files:
-
- 5 edited
-
zdfddm.F90 (modified) (1 diff)
-
zdfevd.F90 (modified) (2 diffs)
-
zdfmxl.F90 (modified) (9 diffs)
-
zdftke.F90 (modified) (15 diffs)
-
zdftmx.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90
r7960 r9987 177 177 & + 0.15 * zrau(ji,jj) * zmskd2(ji,jj) ) 178 178 ! add to the eddy viscosity coef. previously computed 179 # if defined key_zdftmx_new 180 ! key_zdftmx_new: New internal wave-driven param: use avs value computed by zdftmx 181 avs (ji,jj,jk) = avs(ji,jj,jk) + zavfs + zavds 182 # else 179 183 avs (ji,jj,jk) = avt(ji,jj,jk) + zavfs + zavds 184 # endif 180 185 avt (ji,jj,jk) = avt(ji,jj,jk) + zavft + zavdt 181 186 avm (ji,jj,jk) = avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds ) -
branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90
r7960 r9987 19 19 USE zdf_oce ! ocean vertical physics variables 20 20 USE zdfkpp ! KPP vertical mixing 21 USE trd_oce ! trends: ocean variables 22 USE trdtra ! trends manager: tracers 21 23 USE in_out_manager ! I/O manager 22 24 USE iom ! for iom_put … … 122 124 zavt_evd(:,:,:) = avt(:,:,:) - zavt_evd(:,:,:) ! change in avt due to evd 123 125 CALL iom_put( "avt_evd", zavt_evd ) ! output this change 126 IF( l_trdtra ) CALL trd_tra( kt, 'TRA', jp_tem, jptra_evd, zavt_evd ) 124 127 ! 125 128 IF( nn_timing == 1 ) CALL timing_stop('zdf_evd') -
branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90
r7960 r9987 18 18 USE phycst ! physical constants 19 19 USE iom ! I/O library 20 USE eosbn2 ! for zdf_mxl_zint 20 21 USE lib_mpp ! MPP library 21 22 USE wrk_nemo ! work arrays … … 27 28 28 29 PUBLIC zdf_mxl ! called by step.F90 30 PUBLIC zdf_mxl_alloc ! Used in zdf_tke_init 29 31 30 32 INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: nmln !: number of level in the mixed layer (used by TOP) … … 32 34 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmlp !: mixed layer depth (rho=rho0+zdcrit) [m] 33 35 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmlpt !: mixed layer depth at t-points [m] 36 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: hmld_zint !: vertically-interpolated mixed layer depth [m] 37 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: htc_mld ! Heat content of hmld_zint 38 LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ll_found ! Is T_b to be found by interpolation ? 39 LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ll_belowml ! Flag points below mixed layer when ll_found=F 34 40 35 41 REAL(wp), PUBLIC :: rho_c = 0.01_wp !: density criterion for mixed layer depth 36 42 REAL(wp) :: avt_c = 5.e-4_wp ! Kz criterion for the turbocline depth 43 44 TYPE, PUBLIC :: MXL_ZINT !: Structure for MLD defs 45 INTEGER :: mld_type ! mixed layer type 46 REAL(wp) :: zref ! depth of initial T_ref 47 REAL(wp) :: dT_crit ! Critical temp diff 48 REAL(wp) :: iso_frac ! Fraction of rn_dT_crit used 49 END TYPE MXL_ZINT 37 50 38 51 !! * Substitutions … … 51 64 zdf_mxl_alloc = 0 ! set to zero if no array to be allocated 52 65 IF( .NOT. ALLOCATED( nmln ) ) THEN 53 ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), STAT= zdf_mxl_alloc ) 66 ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), hmld_zint(jpi,jpj), & 67 & htc_mld(jpi,jpj), & 68 & ll_found(jpi,jpj), ll_belowml(jpi,jpj,jpk), STAT= zdf_mxl_alloc ) 54 69 ! 55 70 IF( lk_mpp ) CALL mpp_sum ( zdf_mxl_alloc ) … … 79 94 INTEGER, INTENT(in) :: kt ! ocean time-step index 80 95 ! 81 INTEGER :: ji, jj, jk ! dummy loop indices82 INTEGER :: iikn, iiki, ikt , imkt! local integer83 REAL(wp) :: zN2_c ! local scalar96 INTEGER :: ji, jj, jk ! dummy loop indices 97 INTEGER :: iikn, iiki, ikt ! local integer 98 REAL(wp) :: zN2_c ! local scalar 84 99 INTEGER, POINTER, DIMENSION(:,:) :: imld ! 2D workspace 85 100 !!---------------------------------------------------------------------- … … 89 104 CALL wrk_alloc( jpi,jpj, imld ) 90 105 91 IF( kt == nit000 ) THEN106 IF( kt <= nit000 ) THEN 92 107 IF(lwp) WRITE(numout,*) 93 108 IF(lwp) WRITE(numout,*) 'zdf_mxl : mixed layer depth' … … 116 131 DO jj = 1, jpj 117 132 DO ji = 1, jpi 118 imkt = mikt(ji,jj) 119 IF( avt (ji,jj,jk) < avt_c ) imld(ji,jj) = MAX( imkt, jk ) ! Turbocline 133 IF( avt (ji,jj,jk) < avt_c * wmask(ji,jj,jk) ) imld(ji,jj) = jk ! Turbocline 120 134 END DO 121 135 END DO … … 126 140 iiki = imld(ji,jj) 127 141 iikn = nmln(ji,jj) 128 imkt = mikt(ji,jj) 129 hmld (ji,jj) = ( fsdepw(ji,jj,iiki ) - fsdepw(ji,jj,imkt ) ) * ssmask(ji,jj) ! Turbocline depth 130 hmlp (ji,jj) = ( fsdepw(ji,jj,iikn ) - fsdepw(ji,jj,MAX( imkt,nla10 ) ) ) * ssmask(ji,jj) ! Mixed layer depth 131 hmlpt(ji,jj) = ( fsdept(ji,jj,iikn-1) - fsdepw(ji,jj,imkt ) ) * ssmask(ji,jj) ! depth of the last T-point inside the mixed layer 142 hmld (ji,jj) = fsdepw(ji,jj,iiki ) * ssmask(ji,jj) ! Turbocline depth 143 hmlp (ji,jj) = fsdepw(ji,jj,iikn ) * ssmask(ji,jj) ! Mixed layer depth 144 hmlpt(ji,jj) = fsdept(ji,jj,iikn-1) * ssmask(ji,jj) ! depth of the last T-point inside the mixed layer 132 145 END DO 133 146 END DO 134 IF( .NOT.lk_offline ) THEN ! no need to output in offline mode 135 CALL iom_put( "mldr10_1", hmlp ) ! mixed layer depth 136 CALL iom_put( "mldkz5" , hmld ) ! turbocline depth 147 ! no need to output in offline mode 148 IF( .NOT.lk_offline ) THEN 149 IF ( iom_use("mldr10_1") ) THEN 150 IF( ln_isfcav ) THEN 151 CALL iom_put( "mldr10_1", hmlp - risfdep) ! mixed layer thickness 152 ELSE 153 CALL iom_put( "mldr10_1", hmlp ) ! mixed layer depth 154 END IF 155 END IF 156 IF ( iom_use("mldkz5") ) THEN 157 IF( ln_isfcav ) THEN 158 CALL iom_put( "mldkz5" , hmld - risfdep ) ! turbocline thickness 159 ELSE 160 CALL iom_put( "mldkz5" , hmld ) ! turbocline depth 161 END IF 162 END IF 137 163 ENDIF 138 164 165 ! Vertically-interpolated mixed-layer depth diagnostic 166 CALL zdf_mxl_zint( kt ) 167 139 168 IF(ln_ctl) CALL prt_ctl( tab2d_1=REAL(nmln,wp), clinfo1=' nmln : ', tab2d_2=hmlp, clinfo2=' hmlp : ', ovlap=1 ) 140 169 ! … … 144 173 ! 145 174 END SUBROUTINE zdf_mxl 175 176 SUBROUTINE zdf_mxl_zint_mld( sf ) 177 !!---------------------------------------------------------------------------------- 178 !! *** ROUTINE zdf_mxl_zint_mld *** 179 ! 180 ! Calculate vertically-interpolated mixed layer depth diagnostic. 181 ! 182 ! This routine can calculate the mixed layer depth diagnostic suggested by 183 ! Kara et al, 2000, JGR, 105, 16803, but is more general and can calculate 184 ! vertically-interpolated mixed-layer depth diagnostics with other parameter 185 ! settings set in the namzdf_mldzint namelist. 186 ! 187 ! If mld_type=1 the mixed layer depth is calculated as the depth at which the 188 ! density has increased by an amount equivalent to a temperature difference of 189 ! 0.8C at the surface. 190 ! 191 ! For other values of mld_type the mixed layer is calculated as the depth at 192 ! which the temperature differs by 0.8C from the surface temperature. 193 ! 194 ! David Acreman, Daley Calvert 195 ! 196 !!----------------------------------------------------------------------------------- 197 198 TYPE(MXL_ZINT), INTENT(in) :: sf 199 200 ! Diagnostic criteria 201 INTEGER :: nn_mld_type ! mixed layer type 202 REAL(wp) :: rn_zref ! depth of initial T_ref 203 REAL(wp) :: rn_dT_crit ! Critical temp diff 204 REAL(wp) :: rn_iso_frac ! Fraction of rn_dT_crit used 205 206 ! Local variables 207 REAL(wp), PARAMETER :: zepsilon = 1.e-30 ! local small value 208 INTEGER, POINTER, DIMENSION(:,:) :: ikmt ! number of active tracer levels 209 INTEGER, POINTER, DIMENSION(:,:) :: ik_ref ! index of reference level 210 INTEGER, POINTER, DIMENSION(:,:) :: ik_iso ! index of last uniform temp level 211 REAL, POINTER, DIMENSION(:,:,:) :: zT ! Temperature or density 212 REAL, POINTER, DIMENSION(:,:) :: ppzdep ! depth for use in calculating d(rho) 213 REAL, POINTER, DIMENSION(:,:) :: zT_ref ! reference temperature 214 REAL :: zT_b ! base temperature 215 REAL, POINTER, DIMENSION(:,:,:) :: zdTdz ! gradient of zT 216 REAL, POINTER, DIMENSION(:,:,:) :: zmoddT ! Absolute temperature difference 217 REAL :: zdz ! depth difference 218 REAL :: zdT ! temperature difference 219 REAL, POINTER, DIMENSION(:,:) :: zdelta_T ! difference critereon 220 REAL, POINTER, DIMENSION(:,:) :: zRHO1, zRHO2 ! Densities 221 INTEGER :: ji, jj, jk ! loop counter 222 223 !!------------------------------------------------------------------------------------- 224 ! 225 CALL wrk_alloc( jpi, jpj, ikmt, ik_ref, ik_iso) 226 CALL wrk_alloc( jpi, jpj, ppzdep, zT_ref, zdelta_T, zRHO1, zRHO2 ) 227 CALL wrk_alloc( jpi, jpj, jpk, zT, zdTdz, zmoddT ) 228 229 ! Unpack structure 230 nn_mld_type = sf%mld_type 231 rn_zref = sf%zref 232 rn_dT_crit = sf%dT_crit 233 rn_iso_frac = sf%iso_frac 234 235 ! Set the mixed layer depth criterion at each grid point 236 IF( nn_mld_type == 0 ) THEN 237 zdelta_T(:,:) = rn_dT_crit 238 zT(:,:,:) = rhop(:,:,:) 239 ELSE IF( nn_mld_type == 1 ) THEN 240 ppzdep(:,:)=0.0 241 call eos ( tsn(:,:,1,:), ppzdep(:,:), zRHO1(:,:) ) 242 ! Use zT temporarily as a copy of tsn with rn_dT_crit added to SST 243 ! [assumes number of tracers less than number of vertical levels] 244 zT(:,:,1:jpts)=tsn(:,:,1,1:jpts) 245 zT(:,:,jp_tem)=zT(:,:,1)+rn_dT_crit 246 CALL eos( zT(:,:,1:jpts), ppzdep(:,:), zRHO2(:,:) ) 247 zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0 248 ! RHO from eos (2d version) doesn't calculate north or east halo: 249 CALL lbc_lnk( zdelta_T, 'T', 1. ) 250 zT(:,:,:) = rhop(:,:,:) 251 ELSE 252 zdelta_T(:,:) = rn_dT_crit 253 zT(:,:,:) = tsn(:,:,:,jp_tem) 254 END IF 255 256 ! Calculate the gradient of zT and absolute difference for use later 257 DO jk = 1 ,jpk-2 258 zdTdz(:,:,jk) = ( zT(:,:,jk+1) - zT(:,:,jk) ) / fse3w(:,:,jk+1) 259 zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) ) 260 END DO 261 262 ! Find density/temperature at the reference level (Kara et al use 10m). 263 ! ik_ref is the index of the box centre immediately above or at the reference level 264 ! Find rn_zref in the array of model level depths and find the ref 265 ! density/temperature by linear interpolation. 266 DO jk = jpkm1, 2, -1 267 WHERE ( fsdept(:,:,jk) > rn_zref ) 268 ik_ref(:,:) = jk - 1 269 zT_ref(:,:) = zT(:,:,jk-1) + zdTdz(:,:,jk-1) * ( rn_zref - fsdept(:,:,jk-1) ) 270 END WHERE 271 END DO 272 273 ! If the first grid box centre is below the reference level then use the 274 ! top model level to get zT_ref 275 WHERE ( fsdept(:,:,1) > rn_zref ) 276 zT_ref = zT(:,:,1) 277 ik_ref = 1 278 END WHERE 279 280 ! The number of active tracer levels is 1 less than the number of active w levels 281 ikmt(:,:) = mbathy(:,:) - 1 282 283 ! Initialize / reset 284 ll_found(:,:) = .false. 285 286 IF ( rn_iso_frac - zepsilon > 0. ) THEN 287 ! Search for a uniform density/temperature region where adjacent levels 288 ! differ by less than rn_iso_frac * deltaT. 289 ! ik_iso is the index of the last level in the uniform layer 290 ! ll_found indicates whether the mixed layer depth can be found by interpolation 291 ik_iso(:,:) = ik_ref(:,:) 292 DO jj = 1, nlcj 293 DO ji = 1, nlci 294 !CDIR NOVECTOR 295 DO jk = ik_ref(ji,jj), ikmt(ji,jj)-1 296 IF ( zmoddT(ji,jj,jk) > ( rn_iso_frac * zdelta_T(ji,jj) ) ) THEN 297 ik_iso(ji,jj) = jk 298 ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) ) 299 EXIT 300 END IF 301 END DO 302 END DO 303 END DO 304 305 ! Use linear interpolation to find depth of mixed layer base where possible 306 hmld_zint(:,:) = rn_zref 307 DO jj = 1, jpj 308 DO ji = 1, jpi 309 IF (ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0) THEN 310 zdz = abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) ) 311 hmld_zint(ji,jj) = fsdept(ji,jj,ik_iso(ji,jj)) + zdz 312 END IF 313 END DO 314 END DO 315 END IF 316 317 ! If ll_found = .false. then calculate MLD using difference of zdelta_T 318 ! from the reference density/temperature 319 320 ! Prevent this section from working on land points 321 WHERE ( tmask(:,:,1) /= 1.0 ) 322 ll_found = .true. 323 END WHERE 324 325 DO jk=1, jpk 326 ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= zdelta_T(:,:) 327 END DO 328 329 ! Set default value where interpolation cannot be used (ll_found=false) 330 DO jj = 1, jpj 331 DO ji = 1, jpi 332 IF ( .not. ll_found(ji,jj) ) hmld_zint(ji,jj) = fsdept(ji,jj,ikmt(ji,jj)) 333 END DO 334 END DO 335 336 DO jj = 1, jpj 337 DO ji = 1, jpi 338 !CDIR NOVECTOR 339 DO jk = ik_ref(ji,jj)+1, ikmt(ji,jj) 340 IF ( ll_found(ji,jj) ) EXIT 341 IF ( ll_belowml(ji,jj,jk) ) THEN 342 zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * SIGN(1.0, zdTdz(ji,jj,jk-1) ) 343 zdT = zT_b - zT(ji,jj,jk-1) 344 zdz = zdT / zdTdz(ji,jj,jk-1) 345 hmld_zint(ji,jj) = fsdept(ji,jj,jk-1) + zdz 346 EXIT 347 END IF 348 END DO 349 END DO 350 END DO 351 352 hmld_zint(:,:) = hmld_zint(:,:)*tmask(:,:,1) 353 ! 354 CALL wrk_dealloc( jpi, jpj, ikmt, ik_ref, ik_iso) 355 CALL wrk_dealloc( jpi, jpj, ppzdep, zT_ref, zdelta_T, zRHO1, zRHO2 ) 356 CALL wrk_dealloc( jpi,jpj, jpk, zT, zdTdz, zmoddT ) 357 ! 358 END SUBROUTINE zdf_mxl_zint_mld 359 360 SUBROUTINE zdf_mxl_zint_htc( kt ) 361 !!---------------------------------------------------------------------- 362 !! *** ROUTINE zdf_mxl_zint_htc *** 363 !! 364 !! ** Purpose : 365 !! 366 !! ** Method : 367 !!---------------------------------------------------------------------- 368 369 INTEGER, INTENT(in) :: kt ! ocean time-step index 370 371 INTEGER :: ji, jj, jk 372 INTEGER :: ikmax 373 REAL(wp) :: zc, zcoef 374 ! 375 INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ilevel 376 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zthick_0, zthick 377 378 !!---------------------------------------------------------------------- 379 380 IF( .NOT. ALLOCATED(ilevel) ) THEN 381 ALLOCATE( ilevel(jpi,jpj), zthick_0(jpi,jpj), & 382 & zthick(jpi,jpj), STAT=ji ) 383 IF( lk_mpp ) CALL mpp_sum(ji) 384 IF( ji /= 0 ) CALL ctl_stop( 'STOP', 'zdf_mxl_zint_htc : unable to allocate arrays' ) 385 ENDIF 386 387 ! Find last whole model T level above the MLD 388 ilevel(:,:) = 0 389 zthick_0(:,:) = 0._wp 390 391 DO jk = 1, jpkm1 392 DO jj = 1, jpj 393 DO ji = 1, jpi 394 zthick_0(ji,jj) = zthick_0(ji,jj) + fse3t(ji,jj,jk) 395 IF( zthick_0(ji,jj) < hmld_zint(ji,jj) ) ilevel(ji,jj) = jk 396 END DO 397 END DO 398 WRITE(numout,*) 'zthick_0(jk =',jk,') =',zthick_0(2,2) 399 WRITE(numout,*) 'fsdepw(jk+1 =',jk+1,') =',fsdepw(2,2,jk+1) 400 END DO 401 402 ! Surface boundary condition 403 IF( lk_vvl ) THEN ; zthick(:,:) = 0._wp ; htc_mld(:,:) = 0._wp 404 ELSE ; zthick(:,:) = sshn(:,:) ; htc_mld(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1) 405 ENDIF 406 407 ! Deepest whole T level above the MLD 408 ikmax = MIN( MAXVAL( ilevel(:,:) ), jpkm1 ) 409 410 ! Integration down to last whole model T level 411 DO jk = 1, ikmax 412 DO jj = 1, jpj 413 DO ji = 1, jpi 414 zc = fse3t(ji,jj,jk) * REAL( MIN( MAX( 0, ilevel(ji,jj) - jk + 1 ) , 1 ) ) ! 0 below ilevel 415 zthick(ji,jj) = zthick(ji,jj) + zc 416 htc_mld(ji,jj) = htc_mld(ji,jj) + zc * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk) 417 END DO 418 END DO 419 END DO 420 421 ! Subsequent partial T level 422 zthick(:,:) = hmld_zint(:,:) - zthick(:,:) ! remaining thickness to reach MLD 423 424 DO jj = 1, jpj 425 DO ji = 1, jpi 426 htc_mld(ji,jj) = htc_mld(ji,jj) + tsn(ji,jj,ilevel(ji,jj)+1,jp_tem) & 427 & * MIN( fse3t(ji,jj,ilevel(ji,jj)+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel(ji,jj)+1) 428 END DO 429 END DO 430 431 WRITE(numout,*) 'htc_mld(after) =',htc_mld(2,2) 432 433 ! Convert to heat content 434 zcoef = rau0 * rcp 435 htc_mld(:,:) = zcoef * htc_mld(:,:) 436 437 END SUBROUTINE zdf_mxl_zint_htc 438 439 SUBROUTINE zdf_mxl_zint( kt ) 440 !!---------------------------------------------------------------------- 441 !! *** ROUTINE zdf_mxl_zint *** 442 !! 443 !! ** Purpose : 444 !! 445 !! ** Method : 446 !!---------------------------------------------------------------------- 447 448 INTEGER, INTENT(in) :: kt ! ocean time-step index 449 450 INTEGER :: ios 451 INTEGER :: jn 452 453 INTEGER :: nn_mld_diag = 0 ! number of diagnostics 454 455 CHARACTER(len=1) :: cmld 456 457 TYPE(MXL_ZINT) :: sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5 458 TYPE(MXL_ZINT), SAVE, DIMENSION(5) :: mld_diags 459 460 NAMELIST/namzdf_mldzint/ nn_mld_diag, sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5 461 462 !!---------------------------------------------------------------------- 463 464 IF( kt == nit000 ) THEN 465 REWIND( numnam_ref ) ! Namelist namzdf_mldzint in reference namelist 466 READ ( numnam_ref, namzdf_mldzint, IOSTAT = ios, ERR = 901) 467 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in reference namelist', lwp ) 468 469 REWIND( numnam_cfg ) ! Namelist namzdf_mldzint in configuration namelist 470 READ ( numnam_cfg, namzdf_mldzint, IOSTAT = ios, ERR = 902 ) 471 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in configuration namelist', lwp ) 472 IF(lwm) WRITE ( numond, namzdf_mldzint ) 473 474 IF( nn_mld_diag > 5 ) CALL ctl_stop( 'STOP', 'zdf_mxl_ini: Specify no more than 5 MLD definitions' ) 475 476 mld_diags(1) = sn_mld1 477 mld_diags(2) = sn_mld2 478 mld_diags(3) = sn_mld3 479 mld_diags(4) = sn_mld4 480 mld_diags(5) = sn_mld5 481 482 IF( lwp .AND. (nn_mld_diag > 0) ) THEN 483 WRITE(numout,*) '=============== Vertically-interpolated mixed layer ================' 484 WRITE(numout,*) '(Diagnostic number, nn_mld_type, rn_zref, rn_dT_crit, rn_iso_frac)' 485 DO jn = 1, nn_mld_diag 486 WRITE(numout,*) 'MLD criterion',jn,':' 487 WRITE(numout,*) ' nn_mld_type =', mld_diags(jn)%mld_type 488 WRITE(numout,*) ' rn_zref =' , mld_diags(jn)%zref 489 WRITE(numout,*) ' rn_dT_crit =' , mld_diags(jn)%dT_crit 490 WRITE(numout,*) ' rn_iso_frac =', mld_diags(jn)%iso_frac 491 END DO 492 WRITE(numout,*) '====================================================================' 493 ENDIF 494 ENDIF 495 496 IF( nn_mld_diag > 0 ) THEN 497 DO jn = 1, nn_mld_diag 498 WRITE(cmld,'(I1)') jn 499 IF( iom_use( "mldzint_"//cmld ) .OR. iom_use( "mldhtc_"//cmld ) ) THEN 500 CALL zdf_mxl_zint_mld( mld_diags(jn) ) 501 502 IF( iom_use( "mldzint_"//cmld ) ) THEN 503 CALL iom_put( "mldzint_"//cmld, hmld_zint(:,:) ) 504 ENDIF 505 506 IF( iom_use( "mldhtc_"//cmld ) ) THEN 507 CALL zdf_mxl_zint_htc( kt ) 508 CALL iom_put( "mldhtc_"//cmld , htc_mld(:,:) ) 509 ENDIF 510 ENDIF 511 END DO 512 ENDIF 513 514 END SUBROUTINE zdf_mxl_zint 146 515 147 516 !!====================================================================== -
branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
r9188 r9987 53 53 USE timing ! Timing 54 54 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 55 #if defined key_agrif 56 USE agrif_opa_interp 57 USE agrif_opa_update 58 #endif 59 60 55 61 56 62 IMPLICIT NONE … … 77 83 INTEGER :: nn_htau ! type of tke profile of penetration (=0/1) 78 84 REAL(wp) :: rn_efr ! fraction of TKE surface value which penetrates in the ocean 85 REAL(wp) :: rn_c ! fraction of TKE added within the mixed layer by nn_etau 79 86 LOGICAL :: ln_lc ! Langmuir cells (LC) as a source term of TKE or not 80 87 REAL(wp) :: rn_lc ! coef to compute vertical velocity of Langmuir cells … … 82 89 REAL(wp) :: ri_cri ! critic Richardson number (deduced from rn_ediff and rn_ediss values) 83 90 REAL(wp) :: rmxl_min ! minimum mixing length value (deduced from rn_ediff and rn_emin values) [m] 91 REAL(wp) :: rhtau ! coefficient to relate MLD to htau when nn_htau == 2 84 92 REAL(wp) :: rhftau_add = 1.e-3_wp ! add offset applied to HF part of taum (nn_etau=3) 85 93 REAL(wp) :: rhftau_scl = 1.0_wp ! scale factor applied to HF part of taum (nn_etau=3) 86 94 87 95 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: htau ! depth of tke penetration (nn_htau) 96 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_niw !: TKE budget- near-inertial waves term 97 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: efr ! surface boundary condition for nn_etau = 4 88 98 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dissl ! now mixing lenght of dissipation 89 99 #if defined key_c1d … … 108 118 !!---------------------------------------------------------------------- 109 119 ALLOCATE( & 120 & efr (jpi,jpj) , e_niw(jpi,jpj,jpk) , & 110 121 #if defined key_c1d 111 122 & e_dis(jpi,jpj,jpk) , e_mix(jpi,jpj,jpk) , & … … 184 195 avmv_k(:,:,:) = avmv(:,:,:) 185 196 ! 197 #if defined key_agrif 198 ! Update child grid f => parent grid 199 IF( .NOT.Agrif_Root() ) CALL Agrif_Update_Tke( kt ) ! children only 200 #endif 201 ! 186 202 END SUBROUTINE zdf_tke 187 203 … … 312 328 zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) ) 313 329 ! ! TKE Langmuir circulation source term 314 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1) 330 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( 1._wp - fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc ) / & 331 & zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1) 315 332 END DO 316 333 END DO … … 345 362 DO ji = fs_2, fs_jpim1 ! vector opt. 346 363 zcof = zfact1 * tmask(ji,jj,jk) 364 # if defined key_zdftmx_new 365 ! key_zdftmx_new: New internal wave-driven param: set a minimum value for Kz on TKE (ensure numerical stability) 366 zzd_up = zcof * ( MAX( avm(ji,jj,jk+1) + avm(ji,jj,jk), 2.e-5_wp ) ) & ! upper diagonal 367 & / ( fse3t(ji,jj,jk ) * fse3w(ji,jj,jk ) ) 368 zzd_lw = zcof * ( MAX( avm(ji,jj,jk) + avm(ji,jj,jk-1), 2.e-5_wp ) ) & ! lower diagonal 369 & / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk ) ) 370 # else 347 371 zzd_up = zcof * ( avm (ji,jj,jk+1) + avm (ji,jj,jk ) ) & ! upper diagonal 348 372 & / ( fse3t(ji,jj,jk ) * fse3w(ji,jj,jk ) ) 349 373 zzd_lw = zcof * ( avm (ji,jj,jk ) + avm (ji,jj,jk-1) ) & ! lower diagonal 350 374 & / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk ) ) 375 # endif 351 376 ! ! shear prod. at w-point weightened by mask 352 377 zesh2 = ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) & … … 408 433 END DO 409 434 435 ! ! Save TKE prior to nn_etau addition 436 e_niw(:,:,:) = en(:,:,:) 437 ! 410 438 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 411 439 ! ! TKE due to surface and internal wave breaking 412 440 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 441 IF( nn_htau == 2 ) THEN !* mixed-layer depth dependant length scale 442 DO jj = 2, jpjm1 443 DO ji = fs_2, fs_jpim1 ! vector opt. 444 htau(ji,jj) = rhtau * hmlp(ji,jj) 445 END DO 446 END DO 447 ENDIF 448 #if defined key_iomput 449 ! 450 CALL iom_put( "htau", htau(:,:) ) ! Check htau (even if constant in time) 451 #endif 452 ! 413 453 IF( nn_etau == 1 ) THEN !* penetration below the mixed layer (rn_efr fraction) 414 454 DO jk = 2, jpkm1 … … 445 485 END DO 446 486 END DO 487 ELSEIF( nn_etau == 4 ) THEN !* column integral independant of htau (rn_efr must be scaled up) 488 IF( nn_htau == 2 ) THEN ! efr dependant on time-varying htau 489 DO jj = 2, jpjm1 490 DO ji = fs_2, fs_jpim1 ! vector opt. 491 efr(ji,jj) = rn_efr / ( htau(ji,jj) * ( 1._wp - EXP( -bathy(ji,jj) / htau(ji,jj) ) ) ) 492 END DO 493 END DO 494 ENDIF 495 DO jk = 2, jpkm1 496 DO jj = 2, jpjm1 497 DO ji = fs_2, fs_jpim1 ! vector opt. 498 en(ji,jj,jk) = en(ji,jj,jk) + efr(ji,jj) * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) & 499 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) 500 END DO 501 END DO 502 END DO 447 503 ENDIF 448 504 CALL lbc_lnk( en, 'W', 1. ) ! Lateral boundary conditions (sign unchanged) 505 ! 506 DO jk = 2, jpkm1 ! TKE budget: near-inertial waves term 507 DO jj = 2, jpjm1 508 DO ji = fs_2, fs_jpim1 ! vector opt. 509 e_niw(ji,jj,jk) = en(ji,jj,jk) - e_niw(ji,jj,jk) 510 END DO 511 END DO 512 END DO 513 ! 514 CALL lbc_lnk( e_niw, 'W', 1. ) 449 515 ! 450 516 CALL wrk_dealloc( jpi,jpj, imlc ) ! integer … … 705 771 !!---------------------------------------------------------------------- 706 772 INTEGER :: ji, jj, jk ! dummy loop indices 707 INTEGER :: ios 773 INTEGER :: ios, ierr 708 774 !! 709 775 NAMELIST/namzdf_tke/ rn_ediff, rn_ediss , rn_ebb , rn_emin , & 710 776 & rn_emin0, rn_bshear, nn_mxl , ln_mxl0 , & 711 777 & rn_mxl0 , nn_pdl , ln_lc , rn_lc , & 712 & nn_etau , nn_htau , rn_efr 778 & nn_etau , nn_htau , rn_efr , rn_c 713 779 !!---------------------------------------------------------------------- 714 ! 780 715 781 REWIND( numnam_ref ) ! Namelist namzdf_tke in reference namelist : Turbulent Kinetic Energy 716 782 READ ( numnam_ref, namzdf_tke, IOSTAT = ios, ERR = 901) … … 723 789 ! 724 790 ri_cri = 2._wp / ( 2._wp + rn_ediss / rn_ediff ) ! resulting critical Richardson number 791 # if defined key_zdftmx_new 792 ! key_zdftmx_new: New internal wave-driven param: specified value of rn_emin & rmxl_min are used 793 rn_emin = 1.e-10_wp 794 rmxl_min = 1.e-03_wp 795 IF(lwp) THEN ! Control print 796 WRITE(numout,*) 797 WRITE(numout,*) 'zdf_tke_init : New tidal mixing case: force rn_emin = 1.e-10 and rmxl_min = 1.e-3 ' 798 WRITE(numout,*) '~~~~~~~~~~~~' 799 ENDIF 800 # else 725 801 rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) ) ! resulting minimum length to recover molecular viscosity 802 # endif 726 803 ! 727 804 IF(lwp) THEN !* Control print … … 745 822 WRITE(numout,*) ' flag for computation of exp. tke profile nn_htau = ', nn_htau 746 823 WRITE(numout,*) ' fraction of en which pene. the thermocline rn_efr = ', rn_efr 824 WRITE(numout,*) ' fraction of TKE added within the mixed layer by nn_etau rn_c = ', rn_c 747 825 WRITE(numout,*) 748 826 WRITE(numout,*) ' critical Richardson nb with your parameters ri_cri = ', ri_cri … … 755 833 IF( nn_mxl < 0 .OR. nn_mxl > 3 ) CALL ctl_stop( 'bad flag: nn_mxl is 0, 1 or 2 ' ) 756 834 IF( nn_pdl < 0 .OR. nn_pdl > 1 ) CALL ctl_stop( 'bad flag: nn_pdl is 0 or 1 ' ) 757 IF( nn_htau < 0 .OR. nn_htau > 1 ) CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2' )835 IF( nn_htau < 0 .OR. nn_htau > 5 ) CALL ctl_stop( 'bad flag: nn_htau is 0 to 5 ' ) 758 836 IF( nn_etau == 3 .AND. .NOT. ln_cpl ) CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' ) 759 837 … … 763 841 ENDIF 764 842 765 IF( nn_etau == 2 ) CALL zdf_mxl( nit000 ) ! Initialization of nmln 843 IF( nn_etau == 2 ) THEN 844 ierr = zdf_mxl_alloc() 845 nmln(:,:) = nlb10 ! Initialization of nmln 846 ENDIF 847 848 IF( nn_etau /= 0 .and. nn_htau == 2 ) THEN 849 ierr = zdf_mxl_alloc() 850 nmln(:,:) = nlb10 ! Initialization of nmln 851 ENDIF 766 852 767 853 ! !* depth of penetration of surface tke 768 854 IF( nn_etau /= 0 ) THEN 855 htau(:,:) = 0._wp 769 856 SELECT CASE( nn_htau ) ! Choice of the depth of penetration 770 857 CASE( 0 ) ! constant depth penetration (here 10 meters) … … 772 859 CASE( 1 ) ! F(latitude) : 0.5m to 30m poleward of 40 degrees 773 860 htau(:,:) = MAX( 0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) ) ) 861 CASE( 2 ) ! fraction of depth-integrated TKE within mixed-layer 862 rhtau = -1._wp / LOG( 1._wp - rn_c ) 863 CASE( 3 ) ! F(latitude) : 0.5m to 15m poleward of 20 degrees 864 htau(:,:) = MAX( 0.5_wp, MIN( 15._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) ) ) 865 CASE( 4 ) ! F(latitude) : 0.5m to 10m/30m poleward of 13/40 degrees north/south 866 DO jj = 2, jpjm1 867 DO ji = fs_2, fs_jpim1 ! vector opt. 868 IF( gphit(ji,jj) <= 0._wp ) THEN 869 htau(ji,jj) = MAX( 0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) ) ) 870 ELSE 871 htau(ji,jj) = MAX( 0.5_wp, MIN( 10._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) ) ) 872 ENDIF 873 END DO 874 END DO 875 CASE ( 5 ) ! F(latitude) : 0.5m to 10m poleward of 13 degrees north/south, 876 DO jj = 2, jpjm1 ! 10m to 30m between 30/45 degrees south 877 DO ji = fs_2, fs_jpim1 ! vector opt. 878 IF( gphit(ji,jj) <= -30._wp ) THEN 879 htau(ji,jj) = MAX( 10._wp, MIN( 30._wp, 55._wp* ABS( SIN( rpi/120._wp * ( gphit(ji,jj) + 23._wp ) ) ) ) ) 880 ELSE 881 htau(ji,jj) = MAX( 0.5_wp, MIN( 10._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(ji,jj) ) ) ) ) 882 ENDIF 883 END DO 884 END DO 774 885 END SELECT 886 ! 887 IF( nn_etau == 4 .AND. nn_htau /= 2 ) THEN ! efr dependant on constant htau 888 DO jj = 2, jpjm1 889 DO ji = fs_2, fs_jpim1 ! vector opt. 890 efr(ji,jj) = rn_efr / ( htau(ji,jj) * ( 1._wp - EXP( -bathy(ji,jj) / htau(ji,jj) ) ) ) 891 END DO 892 END DO 893 ENDIF 775 894 ENDIF 776 895 ! !* set vertical eddy coef. to the background value -
branches/UKMO/dev_r5518_obs_oper_update_icethick/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90
r7960 r9987 561 561 END SUBROUTINE zdf_tmx_init 562 562 563 #elif defined key_zdftmx_new 564 !!---------------------------------------------------------------------- 565 !! 'key_zdftmx_new' Internal wave-driven vertical mixing 566 !!---------------------------------------------------------------------- 567 !! zdf_tmx : global momentum & tracer Kz with wave induced Kz 568 !! zdf_tmx_init : global momentum & tracer Kz with wave induced Kz 569 !!---------------------------------------------------------------------- 570 USE oce ! ocean dynamics and tracers variables 571 USE dom_oce ! ocean space and time domain variables 572 USE zdf_oce ! ocean vertical physics variables 573 USE zdfddm ! ocean vertical physics: double diffusive mixing 574 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 575 USE eosbn2 ! ocean equation of state 576 USE phycst ! physical constants 577 USE prtctl ! Print control 578 USE in_out_manager ! I/O manager 579 USE iom ! I/O Manager 580 USE lib_mpp ! MPP library 581 USE wrk_nemo ! work arrays 582 USE timing ! Timing 583 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 584 585 IMPLICIT NONE 586 PRIVATE 587 588 PUBLIC zdf_tmx ! called in step module 589 PUBLIC zdf_tmx_init ! called in nemogcm module 590 PUBLIC zdf_tmx_alloc ! called in nemogcm module 591 592 LOGICAL, PUBLIC, PARAMETER :: lk_zdftmx = .TRUE. !: wave-driven mixing flag 593 594 ! !!* Namelist namzdf_tmx : internal wave-driven mixing * 595 INTEGER :: nn_zpyc ! pycnocline-intensified mixing energy proportional to N (=1) or N^2 (=2) 596 LOGICAL :: ln_mevar ! variable (=T) or constant (=F) mixing efficiency 597 LOGICAL :: ln_tsdiff ! account for differential T/S wave-driven mixing (=T) or not (=F) 598 599 REAL(wp) :: r1_6 = 1._wp / 6._wp 600 601 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ebot_tmx ! power available from high-mode wave breaking (W/m2) 602 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: epyc_tmx ! power available from low-mode, pycnocline-intensified wave breaking (W/m2) 603 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ecri_tmx ! power available from low-mode, critical slope wave breaking (W/m2) 604 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hbot_tmx ! WKB decay scale for high-mode energy dissipation (m) 605 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hcri_tmx ! decay scale for low-mode critical slope dissipation (m) 606 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: emix_tmx ! local energy density available for mixing (W/kg) 607 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: bflx_tmx ! buoyancy flux Kz * N^2 (W/kg) 608 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: pcmap_tmx ! vertically integrated buoyancy flux (W/m2) 609 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zav_ratio ! S/T diffusivity ratio (only for ln_tsdiff=T) 610 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zav_wave ! Internal wave-induced diffusivity 611 612 !! * Substitutions 613 # include "zdfddm_substitute.h90" 614 # include "domzgr_substitute.h90" 615 # include "vectopt_loop_substitute.h90" 616 !!---------------------------------------------------------------------- 617 !! NEMO/OPA 4.0 , NEMO Consortium (2016) 618 !! $Id$ 619 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 620 !!---------------------------------------------------------------------- 621 CONTAINS 622 623 INTEGER FUNCTION zdf_tmx_alloc() 624 !!---------------------------------------------------------------------- 625 !! *** FUNCTION zdf_tmx_alloc *** 626 !!---------------------------------------------------------------------- 627 ALLOCATE( ebot_tmx(jpi,jpj), epyc_tmx(jpi,jpj), ecri_tmx(jpi,jpj) , & 628 & hbot_tmx(jpi,jpj), hcri_tmx(jpi,jpj), emix_tmx(jpi,jpj,jpk), & 629 & bflx_tmx(jpi,jpj,jpk), pcmap_tmx(jpi,jpj), zav_ratio(jpi,jpj,jpk), & 630 & zav_wave(jpi,jpj,jpk), STAT=zdf_tmx_alloc ) 631 ! 632 IF( lk_mpp ) CALL mpp_sum ( zdf_tmx_alloc ) 633 IF( zdf_tmx_alloc /= 0 ) CALL ctl_warn('zdf_tmx_alloc: failed to allocate arrays') 634 END FUNCTION zdf_tmx_alloc 635 636 637 SUBROUTINE zdf_tmx( kt ) 638 !!---------------------------------------------------------------------- 639 !! *** ROUTINE zdf_tmx *** 640 !! 641 !! ** Purpose : add to the vertical mixing coefficients the effect of 642 !! breaking internal waves. 643 !! 644 !! ** Method : - internal wave-driven vertical mixing is given by: 645 !! Kz_wave = min( 100 cm2/s, f( Reb = emix_tmx /( Nu * N^2 ) ) 646 !! where emix_tmx is the 3D space distribution of the wave-breaking 647 !! energy and Nu the molecular kinematic viscosity. 648 !! The function f(Reb) is linear (constant mixing efficiency) 649 !! if the namelist parameter ln_mevar = F and nonlinear if ln_mevar = T. 650 !! 651 !! - Compute emix_tmx, the 3D power density that allows to compute 652 !! Reb and therefrom the wave-induced vertical diffusivity. 653 !! This is divided into three components: 654 !! 1. Bottom-intensified low-mode dissipation at critical slopes 655 !! emix_tmx(z) = ( ecri_tmx / rau0 ) * EXP( -(H-z)/hcri_tmx ) 656 !! / ( 1. - EXP( - H/hcri_tmx ) ) * hcri_tmx 657 !! where hcri_tmx is the characteristic length scale of the bottom 658 !! intensification, ecri_tmx a map of available power, and H the ocean depth. 659 !! 2. Pycnocline-intensified low-mode dissipation 660 !! emix_tmx(z) = ( epyc_tmx / rau0 ) * ( sqrt(rn2(z))^nn_zpyc ) 661 !! / SUM( sqrt(rn2(z))^nn_zpyc * e3w(z) ) 662 !! where epyc_tmx is a map of available power, and nn_zpyc 663 !! is the chosen stratification-dependence of the internal wave 664 !! energy dissipation. 665 !! 3. WKB-height dependent high mode dissipation 666 !! emix_tmx(z) = ( ebot_tmx / rau0 ) * rn2(z) * EXP(-z_wkb(z)/hbot_tmx) 667 !! / SUM( rn2(z) * EXP(-z_wkb(z)/hbot_tmx) * e3w(z) ) 668 !! where hbot_tmx is the characteristic length scale of the WKB bottom 669 !! intensification, ebot_tmx is a map of available power, and z_wkb is the 670 !! WKB-stretched height above bottom defined as 671 !! z_wkb(z) = H * SUM( sqrt(rn2(z'>=z)) * e3w(z'>=z) ) 672 !! / SUM( sqrt(rn2(z')) * e3w(z') ) 673 !! 674 !! - update the model vertical eddy viscosity and diffusivity: 675 !! avt = avt + av_wave 676 !! avm = avm + av_wave 677 !! avmu = avmu + mi(av_wave) 678 !! avmv = avmv + mj(av_wave) 679 !! 680 !! - if namelist parameter ln_tsdiff = T, account for differential mixing: 681 !! avs = avt + av_wave * diffusivity_ratio(Reb) 682 !! 683 !! ** Action : - Define emix_tmx used to compute internal wave-induced mixing 684 !! - avt, avs, avm, avmu, avmv increased by internal wave-driven mixing 685 !! 686 !! References : de Lavergne et al. 2015, JPO; 2016, in prep. 687 !!---------------------------------------------------------------------- 688 INTEGER, INTENT(in) :: kt ! ocean time-step 689 ! 690 INTEGER :: ji, jj, jk ! dummy loop indices 691 REAL(wp) :: ztpc ! scalar workspace 692 REAL(wp), DIMENSION(:,:) , POINTER :: zfact ! Used for vertical structure 693 REAL(wp), DIMENSION(:,:) , POINTER :: zhdep ! Ocean depth 694 REAL(wp), DIMENSION(:,:,:), POINTER :: zwkb ! WKB-stretched height above bottom 695 REAL(wp), DIMENSION(:,:,:), POINTER :: zweight ! Weight for high mode vertical distribution 696 REAL(wp), DIMENSION(:,:,:), POINTER :: znu_t ! Molecular kinematic viscosity (T grid) 697 REAL(wp), DIMENSION(:,:,:), POINTER :: znu_w ! Molecular kinematic viscosity (W grid) 698 REAL(wp), DIMENSION(:,:,:), POINTER :: zReb ! Turbulence intensity parameter 699 !!---------------------------------------------------------------------- 700 ! 701 IF( nn_timing == 1 ) CALL timing_start('zdf_tmx') 702 ! 703 CALL wrk_alloc( jpi,jpj, zfact, zhdep ) 704 CALL wrk_alloc( jpi,jpj,jpk, zwkb, zweight, znu_t, znu_w, zReb ) 705 706 ! ! ----------------------------- ! 707 ! ! Internal wave-driven mixing ! (compute zav_wave) 708 ! ! ----------------------------- ! 709 ! 710 ! !* Critical slope mixing: distribute energy over the time-varying ocean depth, 711 ! using an exponential decay from the seafloor. 712 DO jj = 1, jpj ! part independent of the level 713 DO ji = 1, jpi 714 zhdep(ji,jj) = fsdepw(ji,jj,mbkt(ji,jj)+1) ! depth of the ocean 715 zfact(ji,jj) = rau0 * ( 1._wp - EXP( -zhdep(ji,jj) / hcri_tmx(ji,jj) ) ) 716 IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = ecri_tmx(ji,jj) / zfact(ji,jj) 717 END DO 718 END DO 719 720 DO jk = 2, jpkm1 ! complete with the level-dependent part 721 emix_tmx(:,:,jk) = zfact(:,:) * ( EXP( ( fsde3w(:,:,jk ) - zhdep(:,:) ) / hcri_tmx(:,:) ) & 722 & - EXP( ( fsde3w(:,:,jk-1) - zhdep(:,:) ) / hcri_tmx(:,:) ) ) * wmask(:,:,jk) & 723 & / ( fsde3w(:,:,jk) - fsde3w(:,:,jk-1) ) 724 END DO 725 726 ! !* Pycnocline-intensified mixing: distribute energy over the time-varying 727 ! !* ocean depth as proportional to sqrt(rn2)^nn_zpyc 728 729 SELECT CASE ( nn_zpyc ) 730 731 CASE ( 1 ) ! Dissipation scales as N (recommended) 732 733 zfact(:,:) = 0._wp 734 DO jk = 2, jpkm1 ! part independent of the level 735 zfact(:,:) = zfact(:,:) + fse3w(:,:,jk) * SQRT( MAX( 0._wp, rn2(:,:,jk) ) ) * wmask(:,:,jk) 736 END DO 737 738 DO jj = 1, jpj 739 DO ji = 1, jpi 740 IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = epyc_tmx(ji,jj) / ( rau0 * zfact(ji,jj) ) 741 END DO 742 END DO 743 744 DO jk = 2, jpkm1 ! complete with the level-dependent part 745 emix_tmx(:,:,jk) = emix_tmx(:,:,jk) + zfact(:,:) * SQRT( MAX( 0._wp, rn2(:,:,jk) ) ) * wmask(:,:,jk) 746 END DO 747 748 CASE ( 2 ) ! Dissipation scales as N^2 749 750 zfact(:,:) = 0._wp 751 DO jk = 2, jpkm1 ! part independent of the level 752 zfact(:,:) = zfact(:,:) + fse3w(:,:,jk) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk) 753 END DO 754 755 DO jj= 1, jpj 756 DO ji = 1, jpi 757 IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = epyc_tmx(ji,jj) / ( rau0 * zfact(ji,jj) ) 758 END DO 759 END DO 760 761 DO jk = 2, jpkm1 ! complete with the level-dependent part 762 emix_tmx(:,:,jk) = emix_tmx(:,:,jk) + zfact(:,:) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk) 763 END DO 764 765 END SELECT 766 767 ! !* WKB-height dependent mixing: distribute energy over the time-varying 768 ! !* ocean depth as proportional to rn2 * exp(-z_wkb/rn_hbot) 769 770 zwkb(:,:,:) = 0._wp 771 zfact(:,:) = 0._wp 772 DO jk = 2, jpkm1 773 zfact(:,:) = zfact(:,:) + fse3w(:,:,jk) * SQRT( MAX( 0._wp, rn2(:,:,jk) ) ) * wmask(:,:,jk) 774 zwkb(:,:,jk) = zfact(:,:) 775 END DO 776 777 DO jk = 2, jpkm1 778 DO jj = 1, jpj 779 DO ji = 1, jpi 780 IF( zfact(ji,jj) /= 0 ) zwkb(ji,jj,jk) = zhdep(ji,jj) * ( zfact(ji,jj) - zwkb(ji,jj,jk) ) & 781 & * tmask(ji,jj,jk) / zfact(ji,jj) 782 END DO 783 END DO 784 END DO 785 zwkb(:,:,1) = zhdep(:,:) * tmask(:,:,1) 786 787 zweight(:,:,:) = 0._wp 788 DO jk = 2, jpkm1 789 zweight(:,:,jk) = MAX( 0._wp, rn2(:,:,jk) ) * hbot_tmx(:,:) * wmask(:,:,jk) & 790 & * ( EXP( -zwkb(:,:,jk) / hbot_tmx(:,:) ) - EXP( -zwkb(:,:,jk-1) / hbot_tmx(:,:) ) ) 791 END DO 792 793 zfact(:,:) = 0._wp 794 DO jk = 2, jpkm1 ! part independent of the level 795 zfact(:,:) = zfact(:,:) + zweight(:,:,jk) 796 END DO 797 798 DO jj = 1, jpj 799 DO ji = 1, jpi 800 IF( zfact(ji,jj) /= 0 ) zfact(ji,jj) = ebot_tmx(ji,jj) / ( rau0 * zfact(ji,jj) ) 801 END DO 802 END DO 803 804 DO jk = 2, jpkm1 ! complete with the level-dependent part 805 emix_tmx(:,:,jk) = emix_tmx(:,:,jk) + zweight(:,:,jk) * zfact(:,:) * wmask(:,:,jk) & 806 & / ( fsde3w(:,:,jk) - fsde3w(:,:,jk-1) ) 807 END DO 808 809 810 ! Calculate molecular kinematic viscosity 811 znu_t(:,:,:) = 1.e-4_wp * ( 17.91_wp - 0.53810_wp * tsn(:,:,:,jp_tem) + 0.00694_wp * tsn(:,:,:,jp_tem) * tsn(:,:,:,jp_tem) & 812 & + 0.02305_wp * tsn(:,:,:,jp_sal) ) * tmask(:,:,:) * r1_rau0 813 DO jk = 2, jpkm1 814 znu_w(:,:,jk) = 0.5_wp * ( znu_t(:,:,jk-1) + znu_t(:,:,jk) ) * wmask(:,:,jk) 815 END DO 816 817 ! Calculate turbulence intensity parameter Reb 818 DO jk = 2, jpkm1 819 zReb(:,:,jk) = emix_tmx(:,:,jk) / MAX( 1.e-20_wp, znu_w(:,:,jk) * rn2(:,:,jk) ) 820 END DO 821 822 ! Define internal wave-induced diffusivity 823 DO jk = 2, jpkm1 824 zav_wave(:,:,jk) = znu_w(:,:,jk) * zReb(:,:,jk) * r1_6 ! This corresponds to a constant mixing efficiency of 1/6 825 END DO 826 827 IF( ln_mevar ) THEN ! Variable mixing efficiency case : modify zav_wave in the 828 DO jk = 2, jpkm1 ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224 ) regimes 829 DO jj = 1, jpj 830 DO ji = 1, jpi 831 IF( zReb(ji,jj,jk) > 480.00_wp ) THEN 832 zav_wave(ji,jj,jk) = 3.6515_wp * znu_w(ji,jj,jk) * SQRT( zReb(ji,jj,jk) ) 833 ELSEIF( zReb(ji,jj,jk) < 10.224_wp ) THEN 834 zav_wave(ji,jj,jk) = 0.052125_wp * znu_w(ji,jj,jk) * zReb(ji,jj,jk) * SQRT( zReb(ji,jj,jk) ) 835 ENDIF 836 END DO 837 END DO 838 END DO 839 ENDIF 840 841 DO jk = 2, jpkm1 ! Bound diffusivity by molecular value and 100 cm2/s 842 zav_wave(:,:,jk) = MIN( MAX( 1.4e-7_wp, zav_wave(:,:,jk) ), 1.e-2_wp ) * wmask(:,:,jk) 843 END DO 844 845 IF( kt == nit000 ) THEN !* Control print at first time-step: diagnose the energy consumed by zav_wave 846 ztpc = 0._wp 847 DO jk = 2, jpkm1 848 DO jj = 1, jpj 849 DO ji = 1, jpi 850 ztpc = ztpc + fse3w(ji,jj,jk) * e1e2t(ji,jj) & 851 & * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj) 852 END DO 853 END DO 854 END DO 855 IF( lk_mpp ) CALL mpp_sum( ztpc ) 856 ztpc = rau0 * ztpc ! Global integral of rauo * Kz * N^2 = power contributing to mixing 857 858 IF(lwp) THEN 859 WRITE(numout,*) 860 WRITE(numout,*) 'zdf_tmx : Internal wave-driven mixing (tmx)' 861 WRITE(numout,*) '~~~~~~~ ' 862 WRITE(numout,*) 863 WRITE(numout,*) ' Total power consumption by av_wave: ztpc = ', ztpc * 1.e-12_wp, 'TW' 864 ENDIF 865 ENDIF 866 867 ! ! ----------------------- ! 868 ! ! Update mixing coefs ! 869 ! ! ----------------------- ! 870 ! 871 IF( ln_tsdiff ) THEN !* Option for differential mixing of salinity and temperature 872 DO jk = 2, jpkm1 ! Calculate S/T diffusivity ratio as a function of Reb 873 DO jj = 1, jpj 874 DO ji = 1, jpi 875 zav_ratio(ji,jj,jk) = ( 0.505_wp + 0.495_wp * & 876 & TANH( 0.92_wp * ( LOG10( MAX( 1.e-20_wp, zReb(ji,jj,jk) * 5._wp * r1_6 ) ) - 0.60_wp ) ) & 877 & ) * wmask(ji,jj,jk) 878 END DO 879 END DO 880 END DO 881 CALL iom_put( "av_ratio", zav_ratio ) 882 DO jk = 2, jpkm1 !* update momentum & tracer diffusivity with wave-driven mixing 883 fsavs(:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk) * zav_ratio(:,:,jk) 884 avt (:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk) 885 avm (:,:,jk) = avm(:,:,jk) + zav_wave(:,:,jk) 886 END DO 887 ! 888 ELSE !* update momentum & tracer diffusivity with wave-driven mixing 889 DO jk = 2, jpkm1 890 fsavs(:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk) 891 avt (:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk) 892 avm (:,:,jk) = avm(:,:,jk) + zav_wave(:,:,jk) 893 END DO 894 ENDIF 895 896 DO jk = 2, jpkm1 !* update momentum diffusivity at wu and wv points 897 DO jj = 2, jpjm1 898 DO ji = fs_2, fs_jpim1 ! vector opt. 899 avmu(ji,jj,jk) = avmu(ji,jj,jk) + 0.5_wp * ( zav_wave(ji,jj,jk) + zav_wave(ji+1,jj ,jk) ) * wumask(ji,jj,jk) 900 avmv(ji,jj,jk) = avmv(ji,jj,jk) + 0.5_wp * ( zav_wave(ji,jj,jk) + zav_wave(ji ,jj+1,jk) ) * wvmask(ji,jj,jk) 901 END DO 902 END DO 903 END DO 904 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! lateral boundary condition 905 906 ! !* output internal wave-driven mixing coefficient 907 CALL iom_put( "av_wave", zav_wave ) 908 !* output useful diagnostics: N^2, Kz * N^2 (bflx_tmx), 909 ! vertical integral of rau0 * Kz * N^2 (pcmap_tmx), energy density (emix_tmx) 910 IF( iom_use("bflx_tmx") .OR. iom_use("pcmap_tmx") ) THEN 911 bflx_tmx(:,:,:) = MAX( 0._wp, rn2(:,:,:) ) * zav_wave(:,:,:) 912 pcmap_tmx(:,:) = 0._wp 913 DO jk = 2, jpkm1 914 pcmap_tmx(:,:) = pcmap_tmx(:,:) + fse3w(:,:,jk) * bflx_tmx(:,:,jk) * wmask(:,:,jk) 915 END DO 916 pcmap_tmx(:,:) = rau0 * pcmap_tmx(:,:) 917 CALL iom_put( "bflx_tmx", bflx_tmx ) 918 CALL iom_put( "pcmap_tmx", pcmap_tmx ) 919 ENDIF 920 CALL iom_put( "emix_tmx", emix_tmx ) 921 922 CALL wrk_dealloc( jpi,jpj, zfact, zhdep ) 923 CALL wrk_dealloc( jpi,jpj,jpk, zwkb, zweight, znu_t, znu_w, zReb ) 924 925 IF(ln_ctl) CALL prt_ctl(tab3d_1=zav_wave , clinfo1=' tmx - av_wave: ', tab3d_2=avt, clinfo2=' avt: ', ovlap=1, kdim=jpk) 926 ! 927 IF( nn_timing == 1 ) CALL timing_stop('zdf_tmx') 928 ! 929 END SUBROUTINE zdf_tmx 930 931 932 SUBROUTINE zdf_tmx_init 933 !!---------------------------------------------------------------------- 934 !! *** ROUTINE zdf_tmx_init *** 935 !! 936 !! ** Purpose : Initialization of the wave-driven vertical mixing, reading 937 !! of input power maps and decay length scales in netcdf files. 938 !! 939 !! ** Method : - Read the namzdf_tmx namelist and check the parameters 940 !! 941 !! - Read the input data in NetCDF files : 942 !! power available from high-mode wave breaking (mixing_power_bot.nc) 943 !! power available from pycnocline-intensified wave-breaking (mixing_power_pyc.nc) 944 !! power available from critical slope wave-breaking (mixing_power_cri.nc) 945 !! WKB decay scale for high-mode wave-breaking (decay_scale_bot.nc) 946 !! decay scale for critical slope wave-breaking (decay_scale_cri.nc) 947 !! 948 !! ** input : - Namlist namzdf_tmx 949 !! - NetCDF files : mixing_power_bot.nc, mixing_power_pyc.nc, mixing_power_cri.nc, 950 !! decay_scale_bot.nc decay_scale_cri.nc 951 !! 952 !! ** Action : - Increase by 1 the nstop flag is setting problem encounter 953 !! - Define ebot_tmx, epyc_tmx, ecri_tmx, hbot_tmx, hcri_tmx 954 !! 955 !! References : de Lavergne et al. 2015, JPO; 2016, in prep. 956 !! 957 !!---------------------------------------------------------------------- 958 INTEGER :: ji, jj, jk ! dummy loop indices 959 INTEGER :: inum ! local integer 960 INTEGER :: ios 961 REAL(wp) :: zbot, zpyc, zcri ! local scalars 962 !! 963 NAMELIST/namzdf_tmx_new/ nn_zpyc, ln_mevar, ln_tsdiff 964 !!---------------------------------------------------------------------- 965 ! 966 IF( nn_timing == 1 ) CALL timing_start('zdf_tmx_init') 967 ! 968 REWIND( numnam_ref ) ! Namelist namzdf_tmx in reference namelist : Wave-driven mixing 969 READ ( numnam_ref, namzdf_tmx_new, IOSTAT = ios, ERR = 901) 970 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tmx in reference namelist', lwp ) 971 ! 972 REWIND( numnam_cfg ) ! Namelist namzdf_tmx in configuration namelist : Wave-driven mixing 973 READ ( numnam_cfg, namzdf_tmx_new, IOSTAT = ios, ERR = 902 ) 974 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tmx in configuration namelist', lwp ) 975 IF(lwm) WRITE ( numond, namzdf_tmx_new ) 976 ! 977 IF(lwp) THEN ! Control print 978 WRITE(numout,*) 979 WRITE(numout,*) 'zdf_tmx_init : internal wave-driven mixing' 980 WRITE(numout,*) '~~~~~~~~~~~~' 981 WRITE(numout,*) ' Namelist namzdf_tmx_new : set wave-driven mixing parameters' 982 WRITE(numout,*) ' Pycnocline-intensified diss. scales as N (=1) or N^2 (=2) = ', nn_zpyc 983 WRITE(numout,*) ' Variable (T) or constant (F) mixing efficiency = ', ln_mevar 984 WRITE(numout,*) ' Differential internal wave-driven mixing (T) or not (F) = ', ln_tsdiff 985 ENDIF 986 987 ! The new wave-driven mixing parameterization elevates avt and avm in the interior, and 988 ! ensures that avt remains larger than its molecular value (=1.4e-7). Therefore, avtb should 989 ! be set here to a very small value, and avmb to its (uniform) molecular value (=1.4e-6). 990 avmb(:) = 1.4e-6_wp ! viscous molecular value 991 avtb(:) = 1.e-10_wp ! very small diffusive minimum (background avt is specified in zdf_tmx) 992 avtb_2d(:,:) = 1.e0_wp ! uniform 993 IF(lwp) THEN ! Control print 994 WRITE(numout,*) 995 WRITE(numout,*) ' Force the background value applied to avm & avt in TKE to be everywhere ', & 996 & 'the viscous molecular value & a very small diffusive value, resp.' 997 ENDIF 998 999 IF( .NOT.lk_zdfddm ) CALL ctl_stop( 'STOP', 'zdf_tmx_init_new : key_zdftmx_new requires key_zdfddm' ) 1000 1001 ! ! allocate tmx arrays 1002 IF( zdf_tmx_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_tmx_init : unable to allocate tmx arrays' ) 1003 ! 1004 ! ! read necessary fields 1005 CALL iom_open('mixing_power_bot',inum) ! energy flux for high-mode wave breaking [W/m2] 1006 CALL iom_get (inum, jpdom_data, 'field', ebot_tmx, 1 ) 1007 CALL iom_close(inum) 1008 ! 1009 CALL iom_open('mixing_power_pyc',inum) ! energy flux for pynocline-intensified wave breaking [W/m2] 1010 CALL iom_get (inum, jpdom_data, 'field', epyc_tmx, 1 ) 1011 CALL iom_close(inum) 1012 ! 1013 CALL iom_open('mixing_power_cri',inum) ! energy flux for critical slope wave breaking [W/m2] 1014 CALL iom_get (inum, jpdom_data, 'field', ecri_tmx, 1 ) 1015 CALL iom_close(inum) 1016 ! 1017 CALL iom_open('decay_scale_bot',inum) ! spatially variable decay scale for high-mode wave breaking [m] 1018 CALL iom_get (inum, jpdom_data, 'field', hbot_tmx, 1 ) 1019 CALL iom_close(inum) 1020 ! 1021 CALL iom_open('decay_scale_cri',inum) ! spatially variable decay scale for critical slope wave breaking [m] 1022 CALL iom_get (inum, jpdom_data, 'field', hcri_tmx, 1 ) 1023 CALL iom_close(inum) 1024 1025 ebot_tmx(:,:) = ebot_tmx(:,:) * ssmask(:,:) 1026 epyc_tmx(:,:) = epyc_tmx(:,:) * ssmask(:,:) 1027 ecri_tmx(:,:) = ecri_tmx(:,:) * ssmask(:,:) 1028 1029 ! Set once for all to zero the first and last vertical levels of appropriate variables 1030 emix_tmx (:,:, 1 ) = 0._wp 1031 emix_tmx (:,:,jpk) = 0._wp 1032 zav_ratio(:,:, 1 ) = 0._wp 1033 zav_ratio(:,:,jpk) = 0._wp 1034 zav_wave (:,:, 1 ) = 0._wp 1035 zav_wave (:,:,jpk) = 0._wp 1036 1037 zbot = glob_sum( e1e2t(:,:) * ebot_tmx(:,:) ) 1038 zpyc = glob_sum( e1e2t(:,:) * epyc_tmx(:,:) ) 1039 zcri = glob_sum( e1e2t(:,:) * ecri_tmx(:,:) ) 1040 IF(lwp) THEN 1041 WRITE(numout,*) ' High-mode wave-breaking energy: ', zbot * 1.e-12_wp, 'TW' 1042 WRITE(numout,*) ' Pycnocline-intensifed wave-breaking energy: ', zpyc * 1.e-12_wp, 'TW' 1043 WRITE(numout,*) ' Critical slope wave-breaking energy: ', zcri * 1.e-12_wp, 'TW' 1044 ENDIF 1045 ! 1046 IF( nn_timing == 1 ) CALL timing_stop('zdf_tmx_init') 1047 ! 1048 END SUBROUTINE zdf_tmx_init 1049 563 1050 #else 564 1051 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
