- Timestamp:
- 2017-01-20T16:53:40+01:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5518_v3.4_asm_nemovar_community/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90
r6661 r7591 22 22 USE timing ! Timing 23 23 USE trc_oce, ONLY : lk_offline ! offline flag 24 USE lbclnk 25 USE eosbn2 ! Equation of state 24 26 25 27 IMPLICIT NONE … … 35 37 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmlpt !: mixed layer depth at t-points [m] 36 38 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmld_tref !: mixed layer depth at t-points - temperature criterion [m] 39 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hmld_kara !: mixed layer depth of Kara et al [m] 37 40 38 41 REAL(wp), PUBLIC :: rho_c = 0.01_wp !: density criterion for mixed layer depth 39 42 REAL(wp) :: avt_c = 5.e-4_wp ! Kz criterion for the turbocline depth 43 44 ! Namelist variables for namzdf_karaml 45 46 LOGICAL :: ln_kara ! Logical switch for calculating kara mixed 47 ! layer 48 LOGICAL :: ln_kara_write25h ! Logical switch for 25 hour outputs 49 INTEGER :: jpmld_type ! mixed layer type 50 REAL(wp) :: ppz_ref ! depth of initial T_ref 51 REAL(wp) :: ppdT_crit ! Critical temp diff 52 REAL(wp) :: ppiso_frac ! Fraction of ppdT_crit used 53 54 !Used for 25h mean 55 LOGICAL, PRIVATE :: kara_25h_init = .TRUE. !Logical used to initalise 25h 56 !outputs. Necissary, because we need to 57 !initalise the kara_25h on the zeroth 58 !timestep (i.e in the nemogcm_init call) 59 REAL, PRIVATE, ALLOCATABLE, DIMENSION(:,:) :: hmld_kara_25h 40 60 41 61 !! * Substitutions … … 126 146 END DO 127 147 ! depth of the mixing and mixed layers 148 149 CALL zdf_mxl_kara( kt ) ! kara MLD 150 128 151 DO jj = 1, jpj 129 152 DO ji = 1, jpi … … 193 216 END SUBROUTINE zdf_mxl_tref 194 217 218 SUBROUTINE zdf_mxl_kara( kt ) 219 !!---------------------------------------------------------------------------------- 220 !! *** ROUTINE zdf_mxl_kara *** 221 ! 222 ! Calculate mixed layer depth according to the definition of 223 ! Kara et al, 2000, JGR, 105, 16803. 224 ! 225 ! If mld_type=1 the mixed layer depth is calculated as the depth at which the 226 ! density has increased by an amount equivalent to a temperature difference of 227 ! 0.8C at the surface. 228 ! 229 ! For other values of mld_type the mixed layer is calculated as the depth at 230 ! which the temperature differs by 0.8C from the surface temperature. 231 ! 232 ! Original version: David Acreman 233 ! 234 !!----------------------------------------------------------------------------------- 235 236 INTEGER, INTENT( in ) :: kt ! ocean time-step index 237 238 NAMELIST/namzdf_karaml/ ln_kara,jpmld_type, ppz_ref, ppdT_crit, & 239 & ppiso_frac, ln_kara_write25h 240 241 ! Local variables 242 REAL, DIMENSION(jpi,jpk) :: ppzdep ! depth for use in calculating d(rho) 243 REAL(wp), DIMENSION(jpi,jpj,jpts) :: ztsn_2d !Local version of tsn 244 LOGICAL :: ll_found(jpi,jpj) ! Is T_b to be found by interpolation ? 245 LOGICAL :: ll_belowml(jpi,jpj,jpk) ! Flag points below mixed layer when ll_found=F 246 INTEGER :: ji, jj, jk ! loop counter 247 INTEGER :: ik_ref(jpi,jpj) ! index of reference level 248 INTEGER :: ik_iso(jpi,jpj) ! index of last uniform temp level 249 REAL :: zT(jpi,jpj,jpk) ! Temperature or denisty 250 REAL :: zT_ref(jpi,jpj) ! reference temperature 251 REAL :: zT_b ! base temperature 252 REAL :: zdTdz(jpi,jpj,jpk-2) ! gradient of zT 253 REAL :: zmoddT(jpi,jpj,jpk-2) ! Absolute temperature difference 254 REAL :: zdz ! depth difference 255 REAL :: zdT ! temperature difference 256 REAL :: zdelta_T(jpi,jpj) ! difference critereon 257 REAL :: zRHO1(jpi,jpj), zRHO2(jpi,jpj) ! Densities 258 INTEGER, SAVE :: idt, i_steps 259 INTEGER, SAVE :: i_cnt_25h ! Counter for 25 hour means 260 INTEGER :: ios ! Local integer output status for namelist read 261 262 263 264 !!------------------------------------------------------------------------------------- 265 266 IF( kt == nit000 ) THEN 267 ! Default values 268 ln_kara = .FALSE. 269 ln_kara_write25h = .FALSE. 270 jpmld_type = 1 271 ppz_ref = 10.0 272 ppdT_crit = 0.2 273 ppiso_frac = 0.1 274 ! Read namelist 275 REWIND ( numnam_ref ) 276 READ ( numnam_ref, namzdf_karaml, IOSTAT=ios, ERR= 901 ) 277 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_karaml in reference namelist', lwp ) 278 279 REWIND( numnam_cfg ) ! Namelist nam_diadiaop in configuration namelist 3D hourly diagnostics 280 READ ( numnam_cfg, namzdf_karaml, IOSTAT = ios, ERR = 902 ) 281 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_karaml in configuration namelist', lwp ) 282 IF(lwm) WRITE ( numond, namzdf_karaml ) 283 284 285 WRITE(numout,*) '===== Kara mixed layer =====' 286 WRITE(numout,*) 'ln_kara = ', ln_kara 287 WRITE(numout,*) 'jpmld_type = ', jpmld_type 288 WRITE(numout,*) 'ppz_ref = ', ppz_ref 289 WRITE(numout,*) 'ppdT_crit = ', ppdT_crit 290 WRITE(numout,*) 'ppiso_frac = ', ppiso_frac 291 WRITE(numout,*) 'ln_kara_write25h = ', ln_kara_write25h 292 WRITE(numout,*) '============================' 293 294 IF ( .NOT.ln_kara ) THEN 295 WRITE(numout,*) "ln_kara not set; Kara mixed layer not calculated" 296 ELSE 297 IF (.NOT. ALLOCATED(hmld_kara) ) ALLOCATE( hmld_kara(jpi,jpj) ) 298 IF ( ln_kara_write25h .AND. kara_25h_init ) THEN 299 i_cnt_25h = 0 300 IF (.NOT. ALLOCATED(hmld_kara_25h) ) & 301 & ALLOCATE( hmld_kara_25h(jpi,jpj) ) 302 hmld_kara_25h = 0._wp 303 IF( nacc == 1 ) THEN 304 idt = INT(rdtmin) 305 ELSE 306 idt = INT(rdt) 307 ENDIF 308 IF( MOD( 3600,idt) == 0 ) THEN 309 i_steps = 3600 / idt 310 ELSE 311 CALL ctl_stop('STOP', & 312 & 'zdf_mxl_kara: timestep must give MOD(3600,rdt)'// & 313 & ' = 0 otherwise no hourly values are possible') 314 ENDIF 315 ENDIF 316 ENDIF 317 ENDIF 318 319 IF ( ln_kara ) THEN 320 321 !set critical ln_kara 322 ztsn_2d = tsn(:,:,1,:) 323 ztsn_2d(:,:,jp_tem) = ztsn_2d(:,:,jp_tem) + ppdT_crit 324 325 ! Set the mixed layer depth criterion at each grid point 326 ppzdep = 0._wp 327 IF( jpmld_type == 1 ) THEN 328 CALL eos ( tsn(:,:,1,:), & 329 & ppzdep(:,:), zRHO1(:,:) ) 330 CALL eos ( ztsn_2d(:,:,:), & 331 & ppzdep(:,:), zRHO2(:,:) ) 332 zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0 333 ! RHO from eos (2d version) doesn't calculate north or east halo: 334 CALL lbc_lnk( zdelta_T, 'T', 1. ) 335 zT(:,:,:) = rhop(:,:,:) 336 ELSE 337 zdelta_T(:,:) = ppdT_crit 338 zT(:,:,:) = tsn(:,:,:,jp_tem) 339 ENDIF 340 341 ! Calculate the gradient of zT and absolute difference for use later 342 DO jk = 1 ,jpk-2 343 zdTdz(:,:,jk) = ( zT(:,:,jk+1) - zT(:,:,jk) ) / fse3w(:,:,jk+1) 344 zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) ) 345 END DO 346 347 ! Find density/temperature at the reference level (Kara et al use 10m). 348 ! ik_ref is the index of the box centre immediately above or at the reference level 349 ! Find ppz_ref in the array of model level depths and find the ref 350 ! density/temperature by linear interpolation. 351 ik_ref = -1 352 DO jk = jpkm1, 2, -1 353 WHERE( fsdept(:,:,jk) > ppz_ref ) 354 ik_ref(:,:) = jk - 1 355 zT_ref(:,:) = zT(:,:,jk-1) + & 356 & zdTdz(:,:,jk-1) * ( ppz_ref - fsdept(:,:,jk-1) ) 357 ENDWHERE 358 END DO 359 IF ( ANY( ik_ref < 0 ) .OR. ANY( ik_ref > jpkm1 ) ) THEN 360 CALL ctl_stop( "STOP", & 361 & "zdf_mxl_kara: unable to find reference level for kara ML" ) 362 ELSE 363 ! If the first grid box centre is below the reference level then use the 364 ! top model level to get zT_ref 365 WHERE( fsdept(:,:,1) > ppz_ref ) 366 zT_ref = zT(:,:,1) 367 ik_ref = 1 368 ENDWHERE 369 370 ! Search for a uniform density/temperature region where adjacent levels 371 ! differ by less than ppiso_frac * deltaT. 372 ! ik_iso is the index of the last level in the uniform layer 373 ! ll_found indicates whether the mixed layer depth can be found by interpolation 374 ik_iso(:,:) = ik_ref(:,:) 375 ll_found(:,:) = .false. 376 DO jj = 1, nlcj 377 DO ji = 1, nlci 378 !CDIR NOVECTOR 379 DO jk = ik_ref(ji,jj), mbathy(ji,jj)-1 380 IF( zmoddT(ji,jj,jk) > ( ppiso_frac * zdelta_T(ji,jj) ) ) THEN 381 ik_iso(ji,jj) = jk 382 ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) ) 383 EXIT 384 ENDIF 385 END DO 386 END DO 387 END DO 388 389 ! Use linear interpolation to find depth of mixed layer base where possible 390 hmld_kara(:,:) = ppz_ref 391 DO jj = 1, jpj 392 DO ji = 1, jpi 393 IF( ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0 ) THEN 394 zdz = abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) ) 395 hmld_kara(ji,jj) = fsdept(ji,jj,ik_iso(ji,jj)) + zdz 396 ENDIF 397 END DO 398 END DO 399 400 ! If ll_found = .false. then calculate MLD using difference of zdelta_T 401 ! from the reference density/temperature 402 403 ! Prevent this section from working on land points 404 WHERE( tmask(:,:,1) /= 1.0 ) 405 ll_found = .true. 406 ENDWHERE 407 408 DO jk = 1, jpk 409 ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= & 410 & zdelta_T(:,:) 411 END DO 412 413 ! Set default value where interpolation cannot be used (ll_found=false) 414 DO jj = 1, jpj 415 DO ji = 1, jpi 416 IF( .NOT. ll_found(ji,jj) ) & 417 & hmld_kara(ji,jj) = fsdept(ji,jj,mbathy(ji,jj)) 418 END DO 419 END DO 420 421 DO jj = 1, jpj 422 DO ji = 1, jpi 423 !CDIR NOVECTOR 424 DO jk = ik_ref(ji,jj)+1, mbathy(ji,jj) 425 IF( ll_found(ji,jj) ) EXIT 426 IF( ll_belowml(ji,jj,jk) ) THEN 427 zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * & 428 & SIGN(1.0, zdTdz(ji,jj,jk-1) ) 429 zdT = zT_b - zT(ji,jj,jk-1) 430 zdz = zdT / zdTdz(ji,jj,jk-1) 431 hmld_kara(ji,jj) = fsdept(ji,jj,jk-1) + zdz 432 EXIT 433 ENDIF 434 END DO 435 END DO 436 END DO 437 438 hmld_kara(:,:) = hmld_kara(:,:) * tmask(:,:,1) 439 440 IF( ln_kara_write25h ) THEN 441 !Double IF required as i_steps not defined if ln_kara_write25h = 442 ! FALSE 443 IF ( ( MOD( kt, i_steps ) == 0 ) .OR. kara_25h_init ) THEN 444 hmld_kara_25h = hmld_kara_25h + hmld_kara 445 i_cnt_25h = i_cnt_25h + 1 446 IF ( kara_25h_init ) kara_25h_init = .FALSE. 447 ENDIF 448 ENDIF 449 450 #if defined key_iomput 451 IF( kt /= nit000 ) THEN 452 CALL iom_put( "mldkara" , hmld_kara ) 453 IF( ( MOD( i_cnt_25h, 25) == 0 ) .AND. ln_kara_write25h ) & 454 CALL iom_put( "kara25h" , ( hmld_kara_25h / 25._wp ) ) 455 ENDIF 456 #endif 457 458 ENDIF 459 ENDIF 460 461 END SUBROUTINE zdf_mxl_kara 462 195 463 !!====================================================================== 196 464 END MODULE zdfmxl
Note: See TracChangeset
for help on using the changeset viewer.