Changeset 15383
- Timestamp:
- 2021-10-15T13:56:29+02:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE/ZDF/zdfmxl.F90
r15381 r15383 2 2 !!====================================================================== 3 3 !! *** MODULE zdfmxl *** 4 !! Ocean physics: mixed layer depth 4 !! Ocean physics: mixed layer depth 5 5 !!====================================================================== 6 6 !! History : 1.0 ! 2003-08 (G. Madec) original code 7 7 !! 3.2 ! 2009-07 (S. Masson, G. Madec) IOM + merge of DO-loop 8 !! 3.7 ! 2012-03 (G. Madec) make public the density criteria for trdmxl 9 !! - ! 2014-02 (F. Roquet) mixed layer depth calculated using N2 instead of rhop 8 !! 3.7 ! 2012-03 (G. Madec) make public the density criteria for trdmxl 9 !! - ! 2014-02 (F. Roquet) mixed layer depth calculated using N2 instead of rhop 10 10 !!---------------------------------------------------------------------- 11 11 !! zdf_mxl : Compute the turbocline and mixed layer depths. … … 44 44 45 45 46 47 ! Namelist variables for namzdf_karaml 46 47 ! Namelist variables for namzdf_karaml 48 48 LOGICAL :: ln_kara ! Logical switch for calculating kara mixed 49 49 ! layer 50 50 LOGICAL :: ln_kara_write25h ! Logical switch for 25 hour outputs 51 INTEGER :: jpmld_type ! mixed layer type 52 REAL(wp) :: ppz_ref ! depth of initial T_ref 53 REAL(wp) :: ppdT_crit ! Critical temp diff 54 REAL(wp) :: ppiso_frac ! Fraction of ppdT_crit used 55 51 INTEGER :: jpmld_type ! mixed layer type 52 REAL(wp) :: ppz_ref ! depth of initial T_ref 53 REAL(wp) :: ppdT_crit ! Critical temp diff 54 REAL(wp) :: ppiso_frac ! Fraction of ppdT_crit used 55 56 56 !Used for 25h mean 57 LOGICAL, PRIVATE :: kara_25h_init = .TRUE. !Logical used to initalise 25h 57 LOGICAL, PRIVATE :: kara_25h_init = .TRUE. !Logical used to initalise 25h 58 58 !outputs. Necissary, because we need to 59 59 !initalise the kara_25h on the zeroth … … 69 69 70 70 TYPE, PUBLIC :: MXL_ZINT !: Structure for MLD defs 71 INTEGER :: mld_type ! mixed layer type 71 INTEGER :: mld_type ! mixed layer type 72 72 REAL(wp) :: zref ! depth of initial T_ref 73 73 REAL(wp) :: dT_crit ! Critical temp diff 74 REAL(wp) :: iso_frac ! Fraction of rn_dT_crit 74 REAL(wp) :: iso_frac ! Fraction of rn_dT_crit 75 75 END TYPE MXL_ZINT 76 76 … … 89 89 IF( .NOT. ALLOCATED( nmln ) ) THEN 90 90 ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), hmld_zint(jpi,jpj), & 91 & htc_mld(jpi,jpj), ll_found(jpi,jpj), ll_belowml(jpi,jpj,jpk), STAT= zdf_mxl_alloc ) 91 & htc_mld(jpi,jpj), ll_found(jpi,jpj), ll_belowml(jpi,jpj,jpk), STAT= zdf_mxl_alloc ) 92 92 ! 93 93 CALL mpp_sum ( 'zdfmxl', zdf_mxl_alloc ) … … 101 101 !!---------------------------------------------------------------------- 102 102 !! *** ROUTINE zdfmxl *** 103 !! 103 !! 104 104 !! ** Purpose : Compute the turbocline depth and the mixed layer depth 105 105 !! with density criteria. 106 106 !! 107 !! ** Method : The mixed layer depth is the shallowest W depth with 107 !! ** Method : The mixed layer depth is the shallowest W depth with 108 108 !! the density of the corresponding T point (just bellow) bellow a 109 109 !! given value defined locally as rho(10m) + rho_c … … 136 136 zN2_c = grav * rho_c * r1_rau0 ! convert density criteria into N^2 criteria 137 137 DO jk = nlb10, jpkm1 138 DO jj = 1, jpj ! Mixed layer level: w-level 138 DO jj = 1, jpj ! Mixed layer level: w-level 139 139 DO ji = 1, jpi 140 140 ikt = mbkt(ji,jj) … … 147 147 ! w-level of the turbocline and mixing layer (iom_use) 148 148 imld(:,:) = mbkt(:,:) + 1 ! Initialization to the number of w ocean point 149 DO jk = jpkm1, nlb10, -1 ! from the bottom to nlb10 149 DO jk = jpkm1, nlb10, -1 ! from the bottom to nlb10 150 150 DO jj = 1, jpj 151 151 DO ji = 1, jpi 152 IF( avt (ji,jj,jk) < avt_c * wmask(ji,jj,jk) ) imld(ji,jj) = jk ! Turbocline 152 IF( avt (ji,jj,jk) < avt_c * wmask(ji,jj,jk) ) imld(ji,jj) = jk ! Turbocline 153 153 END DO 154 154 END DO … … 163 163 iiki = imld(ji,jj) 164 164 iikn = nmln(ji,jj) 165 hmld (ji,jj) = gdepw_n(ji,jj,iiki ) * ssmask(ji,jj) ! Turbocline depth 165 hmld (ji,jj) = gdepw_n(ji,jj,iiki ) * ssmask(ji,jj) ! Turbocline depth 166 166 hmlp (ji,jj) = gdepw_n(ji,jj,iikn ) * ssmask(ji,jj) ! Mixed layer depth 167 167 hmlpt(ji,jj) = gdept_n(ji,jj,iikn-1) * ssmask(ji,jj) ! depth of the last T-point inside the mixed layer … … 189 189 END SUBROUTINE zdf_mxl 190 190 191 SUBROUTINE zdf_mxl_zint_mld( sf ) 192 !!---------------------------------------------------------------------------------- 193 !! *** ROUTINE zdf_mxl_zint_mld *** 194 ! 195 ! Calculate vertically-interpolated mixed layer depth diagnostic. 196 ! 191 SUBROUTINE zdf_mxl_zint_mld( sf ) 192 !!---------------------------------------------------------------------------------- 193 !! *** ROUTINE zdf_mxl_zint_mld *** 194 ! 195 ! Calculate vertically-interpolated mixed layer depth diagnostic. 196 ! 197 197 ! This routine can calculate the mixed layer depth diagnostic suggested by 198 198 ! Kara et al, 2000, JGR, 105, 16803, but is more general and can calculate 199 199 ! vertically-interpolated mixed-layer depth diagnostics with other parameter 200 ! settings set in the namzdf_mldzint namelist. 201 ! 202 ! If mld_type=1 the mixed layer depth is calculated as the depth at which the 203 ! density has increased by an amount equivalent to a temperature difference of 204 ! 0.8C at the surface. 205 ! 206 ! For other values of mld_type the mixed layer is calculated as the depth at 207 ! which the temperature differs by 0.8C from the surface temperature. 208 ! 209 ! David Acreman, Daley Calvert 210 ! 211 !!----------------------------------------------------------------------------------- 200 ! settings set in the namzdf_mldzint namelist. 201 ! 202 ! If mld_type=1 the mixed layer depth is calculated as the depth at which the 203 ! density has increased by an amount equivalent to a temperature difference of 204 ! 0.8C at the surface. 205 ! 206 ! For other values of mld_type the mixed layer is calculated as the depth at 207 ! which the temperature differs by 0.8C from the surface temperature. 208 ! 209 ! David Acreman, Daley Calvert 210 ! 211 !!----------------------------------------------------------------------------------- 212 212 213 213 TYPE(MXL_ZINT), INTENT(in) :: sf 214 214 215 215 ! Diagnostic criteria 216 INTEGER :: nn_mld_type ! mixed layer type 216 INTEGER :: nn_mld_type ! mixed layer type 217 217 REAL(wp) :: rn_zref ! depth of initial T_ref 218 218 REAL(wp) :: rn_dT_crit ! Critical temp diff … … 221 221 ! Local variables 222 222 REAL(wp), PARAMETER :: zepsilon = 1.e-30 ! local small value 223 INTEGER, DIMENSION(jpi,jpj) :: ikmt ! number of active tracer levels 224 INTEGER, DIMENSION(jpi,jpj) :: ik_ref ! index of reference level 225 INTEGER, DIMENSION(jpi,jpj) :: ik_iso ! index of last uniform temp level 226 REAL, DIMENSION(jpi,jpj,jpk) :: zT ! Temperature or density 227 REAL, DIMENSION(jpi,jpj) :: ppzdep ! depth for use in calculating d(rho) 228 REAL, DIMENSION(jpi,jpj) :: zT_ref ! reference temperature 229 REAL :: zT_b ! base temperature 230 REAL, DIMENSION(jpi,jpj,jpk) :: zdTdz ! gradient of zT 231 REAL, DIMENSION(jpi,jpj,jpk) :: zmoddT ! Absolute temperature difference 232 REAL :: zdz ! depth difference 233 REAL :: zdT ! temperature difference 234 REAL, DIMENSION(jpi,jpj) :: zdelta_T ! difference critereon 235 REAL, DIMENSION(jpi,jpj) :: zRHO1, zRHO2 ! Densities 236 INTEGER :: ji, jj, jk ! loop counter 237 238 !!------------------------------------------------------------------------------------- 239 ! 223 INTEGER, DIMENSION(jpi,jpj) :: ikmt ! number of active tracer levels 224 INTEGER, DIMENSION(jpi,jpj) :: ik_ref ! index of reference level 225 INTEGER, DIMENSION(jpi,jpj) :: ik_iso ! index of last uniform temp level 226 REAL, DIMENSION(jpi,jpj,jpk) :: zT ! Temperature or density 227 REAL, DIMENSION(jpi,jpj) :: ppzdep ! depth for use in calculating d(rho) 228 REAL, DIMENSION(jpi,jpj) :: zT_ref ! reference temperature 229 REAL :: zT_b ! base temperature 230 REAL, DIMENSION(jpi,jpj,jpk) :: zdTdz ! gradient of zT 231 REAL, DIMENSION(jpi,jpj,jpk) :: zmoddT ! Absolute temperature difference 232 REAL :: zdz ! depth difference 233 REAL :: zdT ! temperature difference 234 REAL, DIMENSION(jpi,jpj) :: zdelta_T ! difference critereon 235 REAL, DIMENSION(jpi,jpj) :: zRHO1, zRHO2 ! Densities 236 INTEGER :: ji, jj, jk ! loop counter 237 238 !!------------------------------------------------------------------------------------- 239 ! 240 240 ! Unpack structure 241 241 nn_mld_type = sf%mld_type … … 244 244 rn_iso_frac = sf%iso_frac 245 245 246 ! Set the mixed layer depth criterion at each grid point 246 ! Set the mixed layer depth criterion at each grid point 247 247 IF( nn_mld_type == 0 ) THEN 248 248 zdelta_T(:,:) = rn_dT_crit 249 249 zT(:,:,:) = rhop(:,:,:) 250 250 ELSE IF( nn_mld_type == 1 ) THEN 251 ppzdep(:,:)=0.0 252 call eos ( tsn(:,:,1,:), ppzdep(:,:), zRHO1(:,:) ) 253 ! Use zT temporarily as a copy of tsn with rn_dT_crit added to SST 254 ! [assumes number of tracers less than number of vertical levels] 255 zT(:,:,1:jpts)=tsn(:,:,1,1:jpts) 256 zT(:,:,jp_tem)=zT(:,:,1)+rn_dT_crit 257 CALL eos( zT(:,:,1:jpts), ppzdep(:,:), zRHO2(:,:) ) 258 zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0 259 ! RHO from eos (2d version) doesn't calculate north or east halo: 260 CALL lbc_lnk( 'zdfmxl', zdelta_T, 'T', 1. ) 261 zT(:,:,:) = rhop(:,:,:) 262 ELSE 263 zdelta_T(:,:) = rn_dT_crit 264 zT(:,:,:) = tsn(:,:,:,jp_tem) 265 END IF 266 267 ! Calculate the gradient of zT and absolute difference for use later 268 DO jk = 1 ,jpk-2 269 zdTdz(:,:,jk) = ( zT(:,:,jk+1) - zT(:,:,jk) ) / e3w_n(:,:,jk+1) 270 zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) ) 271 END DO 272 273 ! Find density/temperature at the reference level (Kara et al use 10m). 274 ! ik_ref is the index of the box centre immediately above or at the reference level 275 ! Find rn_zref in the array of model level depths and find the ref 276 ! density/temperature by linear interpolation. 277 DO jk = jpkm1, 2, -1 278 WHERE ( gdept_n(:,:,jk) > rn_zref ) 279 ik_ref(:,:) = jk - 1 280 zT_ref(:,:) = zT(:,:,jk-1) + zdTdz(:,:,jk-1) * ( rn_zref - gdept_n(:,:,jk-1) ) 281 END WHERE 282 END DO 283 284 ! If the first grid box centre is below the reference level then use the 285 ! top model level to get zT_ref 286 WHERE ( gdept_n(:,:,1) > rn_zref ) 287 zT_ref = zT(:,:,1) 288 ik_ref = 1 289 END WHERE 290 291 ! The number of active tracer levels is 1 less than the number of active w levels 292 ikmt(:,:) = mbkt(:,:) - 1 251 ppzdep(:,:)=0.0 252 call eos ( tsn(:,:,1,:), ppzdep(:,:), zRHO1(:,:) ) 253 ! Use zT temporarily as a copy of tsn with rn_dT_crit added to SST 254 ! [assumes number of tracers less than number of vertical levels] 255 zT(:,:,1:jpts)=tsn(:,:,1,1:jpts) 256 zT(:,:,jp_tem)=zT(:,:,1)+rn_dT_crit 257 CALL eos( zT(:,:,1:jpts), ppzdep(:,:), zRHO2(:,:) ) 258 zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0 259 ! RHO from eos (2d version) doesn't calculate north or east halo: 260 CALL lbc_lnk( 'zdfmxl', zdelta_T, 'T', 1. ) 261 zT(:,:,:) = rhop(:,:,:) 262 ELSE 263 zdelta_T(:,:) = rn_dT_crit 264 zT(:,:,:) = tsn(:,:,:,jp_tem) 265 END IF 266 267 ! Calculate the gradient of zT and absolute difference for use later 268 DO jk = 1 ,jpk-2 269 zdTdz(:,:,jk) = ( zT(:,:,jk+1) - zT(:,:,jk) ) / e3w_n(:,:,jk+1) 270 zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) ) 271 END DO 272 273 ! Find density/temperature at the reference level (Kara et al use 10m). 274 ! ik_ref is the index of the box centre immediately above or at the reference level 275 ! Find rn_zref in the array of model level depths and find the ref 276 ! density/temperature by linear interpolation. 277 DO jk = jpkm1, 2, -1 278 WHERE ( gdept_n(:,:,jk) > rn_zref ) 279 ik_ref(:,:) = jk - 1 280 zT_ref(:,:) = zT(:,:,jk-1) + zdTdz(:,:,jk-1) * ( rn_zref - gdept_n(:,:,jk-1) ) 281 END WHERE 282 END DO 283 284 ! If the first grid box centre is below the reference level then use the 285 ! top model level to get zT_ref 286 WHERE ( gdept_n(:,:,1) > rn_zref ) 287 zT_ref = zT(:,:,1) 288 ik_ref = 1 289 END WHERE 290 291 ! The number of active tracer levels is 1 less than the number of active w levels 292 ikmt(:,:) = mbkt(:,:) - 1 293 293 294 294 ! Initialize / reset … … 296 296 297 297 IF ( rn_iso_frac - zepsilon > 0. ) THEN 298 ! Search for a uniform density/temperature region where adjacent levels 299 ! differ by less than rn_iso_frac * deltaT. 300 ! ik_iso is the index of the last level in the uniform layer 301 ! ll_found indicates whether the mixed layer depth can be found by interpolation 302 ik_iso(:,:) = ik_ref(:,:) 303 DO jj = 1, nlcj 304 DO ji = 1, nlci 305 !CDIR NOVECTOR 306 DO jk = ik_ref(ji,jj), ikmt(ji,jj)-1 307 IF ( zmoddT(ji,jj,jk) > ( rn_iso_frac * zdelta_T(ji,jj) ) ) THEN 308 ik_iso(ji,jj) = jk 309 ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) ) 310 EXIT 311 END IF 312 END DO 313 END DO 314 END DO 315 316 ! Use linear interpolation to find depth of mixed layer base where possible 317 hmld_zint(:,:) = rn_zref 318 DO jj = 1, jpj 319 DO ji = 1, jpi 320 IF (ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0) THEN 321 zdz = abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) ) 322 hmld_zint(ji,jj) = gdept_n(ji,jj,ik_iso(ji,jj)) + zdz 323 END IF 324 END DO 325 END DO 298 ! Search for a uniform density/temperature region where adjacent levels 299 ! differ by less than rn_iso_frac * deltaT. 300 ! ik_iso is the index of the last level in the uniform layer 301 ! ll_found indicates whether the mixed layer depth can be found by interpolation 302 ik_iso(:,:) = ik_ref(:,:) 303 DO jj = 1, nlcj 304 DO ji = 1, nlci 305 !CDIR NOVECTOR 306 DO jk = ik_ref(ji,jj), ikmt(ji,jj)-1 307 IF ( zmoddT(ji,jj,jk) > ( rn_iso_frac * zdelta_T(ji,jj) ) ) THEN 308 ik_iso(ji,jj) = jk 309 ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) ) 310 EXIT 311 END IF 312 END DO 313 END DO 314 END DO 315 316 ! Use linear interpolation to find depth of mixed layer base where possible 317 hmld_zint(:,:) = rn_zref 318 DO jj = 1, jpj 319 DO ji = 1, jpi 320 IF (ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0) THEN 321 zdz = abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) ) 322 hmld_zint(ji,jj) = gdept_n(ji,jj,ik_iso(ji,jj)) + zdz 323 END IF 324 END DO 325 END DO 326 326 END IF 327 327 328 ! If ll_found = .false. then calculate MLD using difference of zdelta_T 329 ! from the reference density/temperature 330 331 ! Prevent this section from working on land points 332 WHERE ( tmask(:,:,1) /= 1.0 ) 333 ll_found = .true. 334 END WHERE 335 336 DO jk=1, jpk 337 ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= zdelta_T(:,:) 338 END DO 339 340 ! Set default value where interpolation cannot be used (ll_found=false) 341 DO jj = 1, jpj 342 DO ji = 1, jpi 343 IF ( .not. ll_found(ji,jj) ) hmld_zint(ji,jj) = gdept_n(ji,jj,ikmt(ji,jj)) 344 END DO 345 END DO 346 347 DO jj = 1, jpj 348 DO ji = 1, jpi 349 !CDIR NOVECTOR 350 DO jk = ik_ref(ji,jj)+1, ikmt(ji,jj) 351 IF ( ll_found(ji,jj) ) EXIT 352 IF ( ll_belowml(ji,jj,jk) ) THEN 353 zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * SIGN(1.0, zdTdz(ji,jj,jk-1) ) 354 zdT = zT_b - zT(ji,jj,jk-1) 355 zdz = zdT / zdTdz(ji,jj,jk-1) 356 hmld_zint(ji,jj) = gdept_n(ji,jj,jk-1) + zdz 357 EXIT 358 END IF 359 END DO 360 END DO 361 END DO 362 363 hmld_zint(:,:) = hmld_zint(:,:)*tmask(:,:,1) 364 ! 328 ! If ll_found = .false. then calculate MLD using difference of zdelta_T 329 ! from the reference density/temperature 330 331 ! Prevent this section from working on land points 332 WHERE ( tmask(:,:,1) /= 1.0 ) 333 ll_found = .true. 334 END WHERE 335 336 DO jk=1, jpk 337 ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= zdelta_T(:,:) 338 END DO 339 340 ! Set default value where interpolation cannot be used (ll_found=false) 341 DO jj = 1, jpj 342 DO ji = 1, jpi 343 IF ( .not. ll_found(ji,jj) ) hmld_zint(ji,jj) = gdept_n(ji,jj,ikmt(ji,jj)) 344 END DO 345 END DO 346 347 DO jj = 1, jpj 348 DO ji = 1, jpi 349 !CDIR NOVECTOR 350 DO jk = ik_ref(ji,jj)+1, ikmt(ji,jj) 351 IF ( ll_found(ji,jj) ) EXIT 352 IF ( ll_belowml(ji,jj,jk) ) THEN 353 zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * SIGN(1.0, zdTdz(ji,jj,jk-1) ) 354 zdT = zT_b - zT(ji,jj,jk-1) 355 zdz = zdT / zdTdz(ji,jj,jk-1) 356 hmld_zint(ji,jj) = gdept_n(ji,jj,jk-1) + zdz 357 EXIT 358 END IF 359 END DO 360 END DO 361 END DO 362 363 hmld_zint(:,:) = hmld_zint(:,:)*tmask(:,:,1) 364 ! 365 365 END SUBROUTINE zdf_mxl_zint_mld 366 366 … … 368 368 !!---------------------------------------------------------------------- 369 369 !! *** ROUTINE zdf_mxl_zint_htc *** 370 !!371 !! ** Purpose :372 370 !! 373 !! ** Method : 371 !! ** Purpose : 372 !! 373 !! ** Method : 374 374 !!---------------------------------------------------------------------- 375 375 … … 396 396 zthick_0(:,:) = 0._wp 397 397 398 DO jk = 1, jpkm1 398 DO jk = 1, jpkm1 399 399 DO jj = 1, jpj 400 DO ji = 1, jpi 400 DO ji = 1, jpi 401 401 zthick_0(ji,jj) = zthick_0(ji,jj) + e3t_n(ji,jj,jk) 402 402 IF( zthick_0(ji,jj) < hmld_zint(ji,jj) ) ilevel(ji,jj) = jk … … 408 408 409 409 ! Surface boundary condition 410 IF( ln_linssh ) THEN ; zthick(:,:) = sshn(:,:) ; htc_mld(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1) 411 ELSE ; zthick(:,:) = 0._wp ; htc_mld(:,:) = 0._wp 410 IF( ln_linssh ) THEN ; zthick(:,:) = sshn(:,:) ; htc_mld(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1) 411 ELSE ; zthick(:,:) = 0._wp ; htc_mld(:,:) = 0._wp 412 412 ENDIF 413 413 … … 431 431 DO jj = 1, jpj 432 432 DO ji = 1, jpi 433 htc_mld(ji,jj) = htc_mld(ji,jj) + tsn(ji,jj,ilevel(ji,jj)+1,jp_tem) & 433 htc_mld(ji,jj) = htc_mld(ji,jj) + tsn(ji,jj,ilevel(ji,jj)+1,jp_tem) & 434 434 & * MIN( e3t_n(ji,jj,ilevel(ji,jj)+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel(ji,jj)+1) 435 435 END DO … … 447 447 !!---------------------------------------------------------------------- 448 448 !! *** ROUTINE zdf_mxl_zint *** 449 !!450 !! ** Purpose :451 449 !! 452 !! ** Method : 450 !! ** Purpose : 451 !! 452 !! ** Method : 453 453 !!---------------------------------------------------------------------- 454 454 … … 468 468 469 469 !!---------------------------------------------------------------------- 470 470 471 471 IF( kt == nit000 ) THEN 472 REWIND( numnam_ref ) ! Namelist namzdf_mldzint in reference namelist 472 REWIND( numnam_ref ) ! Namelist namzdf_mldzint in reference namelist 473 473 READ ( numnam_ref, namzdf_mldzint, IOSTAT = ios, ERR = 901) 474 474 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in reference namelist' ) 475 475 476 REWIND( numnam_cfg ) ! Namelist namzdf_mldzint in configuration namelist 476 REWIND( numnam_cfg ) ! Namelist namzdf_mldzint in configuration namelist 477 477 READ ( numnam_cfg, namzdf_mldzint, IOSTAT = ios, ERR = 902 ) 478 478 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in configuration namelist' ) … … 523 523 524 524 525 SUBROUTINE zdf_mxl_kara( kt ) 526 !!---------------------------------------------------------------------------------- 527 !! *** ROUTINE zdf_mxl_kara *** 528 ! 529 ! Calculate mixed layer depth according to the definition of 530 ! Kara et al, 2000, JGR, 105, 16803. 531 ! 532 ! If mld_type=1 the mixed layer depth is calculated as the depth at which the 533 ! density has increased by an amount equivalent to a temperature difference of 534 ! 0.8C at the surface. 535 ! 536 ! For other values of mld_type the mixed layer is calculated as the depth at 537 ! which the temperature differs by 0.8C from the surface temperature. 538 ! 539 ! Original version: David Acreman 540 ! 525 SUBROUTINE zdf_mxl_kara( kt ) 526 !!---------------------------------------------------------------------------------- 527 !! *** ROUTINE zdf_mxl_kara *** 528 ! 529 ! Calculate mixed layer depth according to the definition of 530 ! Kara et al, 2000, JGR, 105, 16803. 531 ! 532 ! If mld_type=1 the mixed layer depth is calculated as the depth at which the 533 ! density has increased by an amount equivalent to a temperature difference of 534 ! 0.8C at the surface. 535 ! 536 ! For other values of mld_type the mixed layer is calculated as the depth at 537 ! which the temperature differs by 0.8C from the surface temperature. 538 ! 539 ! Original version: David Acreman 540 ! 541 541 !!----------------------------------------------------------------------------------- 542 543 INTEGER, INTENT( in ) :: kt ! ocean time-step index 544 542 543 INTEGER, INTENT( in ) :: kt ! ocean time-step index 544 545 545 NAMELIST/namzdf_karaml/ ln_kara,jpmld_type, ppz_ref, ppdT_crit, & 546 & ppiso_frac, ln_kara_write25h 547 548 ! Local variables 549 REAL, DIMENSION(jpi,jpk) :: ppzdep ! depth for use in calculating d(rho) 550 REAL(wp), DIMENSION(jpi,jpj,jpts) :: ztsn_2d !Local version of tsn 551 LOGICAL :: ll_found(jpi,jpj) ! Is T_b to be found by interpolation ? 552 LOGICAL :: ll_belowml(jpi,jpj,jpk) ! Flag points below mixed layer when ll_found=F 553 INTEGER :: ji, jj, jk ! loop counter 554 INTEGER :: ik_ref(jpi,jpj) ! index of reference level 555 INTEGER :: ik_iso(jpi,jpj) ! index of last uniform temp level 556 REAL :: zT(jpi,jpj,jpk) ! Temperature or denisty 557 REAL :: zT_ref(jpi,jpj) ! reference temperature 558 REAL :: zT_b ! base temperature 559 REAL :: zdTdz(jpi,jpj,jpk-2) ! gradient of zT 560 REAL :: zmoddT(jpi,jpj,jpk-2) ! Absolute temperature difference 561 REAL :: zdz ! depth difference 562 REAL :: zdT ! temperature difference 563 REAL :: zdelta_T(jpi,jpj) ! difference critereon 546 & ppiso_frac, ln_kara_write25h 547 548 ! Local variables 549 REAL, DIMENSION(jpi,jpk) :: ppzdep ! depth for use in calculating d(rho) 550 REAL(wp), DIMENSION(jpi,jpj,jpts) :: ztsn_2d !Local version of tsn 551 INTEGER, DIMENSION(jpi,jpj) :: ikmt ! number of active tracer levels 552 LOGICAL :: ll_found(jpi,jpj) ! Is T_b to be found by interpolation ? 553 LOGICAL :: ll_belowml(jpi,jpj,jpk) ! Flag points below mixed layer when ll_found=F 554 INTEGER :: ji, jj, jk ! loop counter 555 INTEGER :: ik_ref(jpi,jpj) ! index of reference level 556 INTEGER :: ik_iso(jpi,jpj) ! index of last uniform temp level 557 REAL :: zT(jpi,jpj,jpk) ! Temperature or denisty 558 REAL :: zT_ref(jpi,jpj) ! reference temperature 559 REAL :: zT_b ! base temperature 560 REAL :: zdTdz(jpi,jpj,jpk-2) ! gradient of zT 561 REAL :: zmoddT(jpi,jpj,jpk-2) ! Absolute temperature difference 562 REAL :: zdz ! depth difference 563 REAL :: zdT ! temperature difference 564 REAL :: zdelta_T(jpi,jpj) ! difference critereon 564 565 REAL :: zRHO1(jpi,jpj), zRHO2(jpi,jpj) ! Densities 565 566 INTEGER, SAVE :: idt, i_steps … … 567 568 INTEGER :: ios ! Local integer output status for namelist read 568 569 569 570 571 !!------------------------------------------------------------------------------------- 572 573 IF( kt == nit000 ) THEN 574 ! Default values 570 571 572 !!------------------------------------------------------------------------------------- 573 574 IF( kt == nit000 ) THEN 575 ! Default values 575 576 ln_kara = .FALSE. 576 ln_kara_write25h = .FALSE. 577 jpmld_type = 1 578 ppz_ref = 10.0 579 ppdT_crit = 0.2 580 ppiso_frac = 0.1 577 ln_kara_write25h = .FALSE. 578 jpmld_type = 1 579 ppz_ref = 10.0 580 ppdT_crit = 0.2 581 ppiso_frac = 0.1 581 582 ! Read namelist 582 583 REWIND ( numnam_ref ) … … 588 589 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namzdf_karaml in configuration namelist' ) 589 590 IF(lwm) WRITE ( numond, namzdf_karaml ) 590 591 592 WRITE(numout,*) '===== Kara mixed layer =====' 591 592 593 WRITE(numout,*) '===== Kara mixed layer =====' 593 594 WRITE(numout,*) 'ln_kara = ', ln_kara 594 WRITE(numout,*) 'jpmld_type = ', jpmld_type 595 WRITE(numout,*) 'ppz_ref = ', ppz_ref 596 WRITE(numout,*) 'ppdT_crit = ', ppdT_crit 595 WRITE(numout,*) 'jpmld_type = ', jpmld_type 596 WRITE(numout,*) 'ppz_ref = ', ppz_ref 597 WRITE(numout,*) 'ppdT_crit = ', ppdT_crit 597 598 WRITE(numout,*) 'ppiso_frac = ', ppiso_frac 598 599 WRITE(numout,*) 'ln_kara_write25h = ', ln_kara_write25h 599 WRITE(numout,*) '============================' 600 600 WRITE(numout,*) '============================' 601 601 602 IF ( .NOT.ln_kara ) THEN 602 WRITE(numout,*) "ln_kara not set; Kara mixed layer not calculated" 603 WRITE(numout,*) "ln_kara not set; Kara mixed layer not calculated" 603 604 ELSE 604 605 IF (.NOT. ALLOCATED(hmld_kara) ) ALLOCATE( hmld_kara(jpi,jpj) ) … … 615 616 616 617 idt = INT(rdt) 617 IF( MOD( 3600,idt) == 0 ) THEN 618 i_steps = 3600 / idt 619 ELSE 618 IF( MOD( 3600,idt) == 0 ) THEN 619 i_steps = 3600 / idt 620 ELSE 620 621 CALL ctl_stop('STOP', & 621 622 & 'zdf_mxl_kara: timestep must give MOD(3600,rdt)'// & 622 & ' = 0 otherwise no hourly values are possible') 623 ENDIF 623 & ' = 0 otherwise no hourly values are possible') 624 ENDIF 624 625 ENDIF 625 626 ENDIF 626 627 ENDIF 627 628 628 629 IF ( ln_kara ) THEN 629 630 630 631 !set critical ln_kara 631 632 ztsn_2d = tsn(:,:,1,:) 632 633 ztsn_2d(:,:,jp_tem) = ztsn_2d(:,:,jp_tem) + ppdT_crit 633 634 ! Set the mixed layer depth criterion at each grid point 634 635 ! Set the mixed layer depth criterion at each grid point 635 636 ppzdep = 0._wp 636 IF( jpmld_type == 1 ) THEN 637 IF( jpmld_type == 1 ) THEN 637 638 CALL eos ( tsn(:,:,1,:), & 638 & ppzdep(:,:), zRHO1(:,:) ) 639 & ppzdep(:,:), zRHO1(:,:) ) 639 640 CALL eos ( ztsn_2d(:,:,:), & 640 & ppzdep(:,:), zRHO2(:,:) ) 641 zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0 642 ! RHO from eos (2d version) doesn't calculate north or east halo: 643 CALL lbc_lnk( 'zdf_mxl_kara',zdelta_T, 'T', 1. ) 644 zT(:,:,:) = rhop(:,:,:) 645 ELSE 646 zdelta_T(:,:) = ppdT_crit 647 zT(:,:,:) = tsn(:,:,:,jp_tem) 648 ENDIF 649 650 ! Calculate the gradient of zT and absolute difference for use later 651 DO jk = 1 ,jpk-2 652 zdTdz(:,:,jk) = ( zT(:,:,jk+1) - zT(:,:,jk) ) / e3w_n(:,:,jk+1) 653 zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) ) 654 END DO 655 656 ! Find density/temperature at the reference level (Kara et al use 10m). 657 ! ik_ref is the index of the box centre immediately above or at the reference level 658 ! Find ppz_ref in the array of model level depths and find the ref 659 ! density/temperature by linear interpolation. 641 & ppzdep(:,:), zRHO2(:,:) ) 642 zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0 643 ! RHO from eos (2d version) doesn't calculate north or east halo: 644 CALL lbc_lnk( 'zdf_mxl_kara',zdelta_T, 'T', 1. ) 645 zT(:,:,:) = rhop(:,:,:) 646 ELSE 647 zdelta_T(:,:) = ppdT_crit 648 zT(:,:,:) = tsn(:,:,:,jp_tem) 649 ENDIF 650 651 ! Calculate the gradient of zT and absolute difference for use later 652 DO jk = 1 ,jpk-2 653 zdTdz(:,:,jk) = ( zT(:,:,jk+1) - zT(:,:,jk) ) / e3w_n(:,:,jk+1) 654 zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) ) 655 END DO 656 657 ! Find density/temperature at the reference level (Kara et al use 10m). 658 ! ik_ref is the index of the box centre immediately above or at the reference level 659 ! Find ppz_ref in the array of model level depths and find the ref 660 ! density/temperature by linear interpolation. 660 661 ik_ref = -1 661 DO jk = jpkm1, 2, -1 662 WHERE( gdept_n(:,:,jk) > ppz_ref ) 663 ik_ref(:,:) = jk - 1 662 DO jk = jpkm1, 2, -1 663 WHERE( gdept_n(:,:,jk) > ppz_ref ) 664 ik_ref(:,:) = jk - 1 664 665 zT_ref(:,:) = zT(:,:,jk-1) + & 665 & zdTdz(:,:,jk-1) * ( ppz_ref - gdept_n(:,:,jk-1) ) 666 ENDWHERE 666 & zdTdz(:,:,jk-1) * ( ppz_ref - gdept_n(:,:,jk-1) ) 667 ENDWHERE 667 668 END DO 668 669 IF ( ANY( ik_ref < 0 ) .OR. ANY( ik_ref > jpkm1 ) ) THEN 669 670 CALL ctl_stop( "STOP", & 670 & "zdf_mxl_kara: unable to find reference level for kara ML" ) 671 & "zdf_mxl_kara: unable to find reference level for kara ML" ) 671 672 ELSE 672 ! If the first grid box centre is below the reference level then use the 673 ! top model level to get zT_ref 674 WHERE( gdept_n(:,:,1) > ppz_ref ) 675 zT_ref = zT(:,:,1) 676 ik_ref = 1 677 ENDWHERE 678 679 ! Search for a uniform density/temperature region where adjacent levels 680 ! differ by less than ppiso_frac * deltaT. 681 ! ik_iso is the index of the last level in the uniform layer 682 ! ll_found indicates whether the mixed layer depth can be found by interpolation 683 ik_iso(:,:) = ik_ref(:,:) 684 ll_found(:,:) = .false. 685 DO jj = 1, nlcj 686 DO ji = 1, nlci 687 !CDIR NOVECTOR 688 DO jk = ik_ref(ji,jj), mbkt(ji,jj)-1 !mbathy(ji,jj)-1 689 IF( zmoddT(ji,jj,jk) > ( ppiso_frac * zdelta_T(ji,jj) ) ) THEN 690 ik_iso(ji,jj) = jk 691 ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) ) 692 EXIT 693 ENDIF 694 END DO 695 END DO 696 END DO 697 698 ! Use linear interpolation to find depth of mixed layer base where possible 699 hmld_kara(:,:) = ppz_ref 700 DO jj = 1, jpj 701 DO ji = 1, jpi 702 IF( ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0 ) THEN 703 zdz = abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) ) 704 hmld_kara(ji,jj) = gdept_n(ji,jj,ik_iso(ji,jj)) + zdz 705 ENDIF 706 END DO 707 END DO 708 709 ! If ll_found = .false. then calculate MLD using difference of zdelta_T 710 ! from the reference density/temperature 711 712 ! Prevent this section from working on land points 713 WHERE( tmask(:,:,1) /= 1.0 ) 714 ll_found = .true. 715 ENDWHERE 716 717 DO jk = 1, jpk 673 ! If the first grid box centre is below the reference level then use the 674 ! top model level to get zT_ref 675 WHERE( gdept_n(:,:,1) > ppz_ref ) 676 zT_ref = zT(:,:,1) 677 ik_ref = 1 678 ENDWHERE 679 680 ! The number of active tracer levels is 1 less than the number of active w levels 681 ikmt(:,:) = mbkt(:,:) - 1 682 683 ! Search for a uniform density/temperature region where adjacent levels 684 ! differ by less than ppiso_frac * deltaT. 685 ! ik_iso is the index of the last level in the uniform layer 686 ! ll_found indicates whether the mixed layer depth can be found by interpolation 687 ik_iso(:,:) = ik_ref(:,:) 688 ll_found(:,:) = .false. 689 DO jj = 1, nlcj 690 DO ji = 1, nlci 691 !CDIR NOVECTOR 692 DO jk = ik_ref(ji,jj), ikmt(ji,jj)-1 ! mbkt(ji,jj)-1 !mbathy(ji,jj)-1 693 IF( zmoddT(ji,jj,jk) > ( ppiso_frac * zdelta_T(ji,jj) ) ) THEN 694 ik_iso(ji,jj) = jk 695 ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) ) 696 EXIT 697 ENDIF 698 END DO 699 END DO 700 END DO 701 702 ! Use linear interpolation to find depth of mixed layer base where possible 703 hmld_kara(:,:) = ppz_ref 704 DO jj = 1, jpj 705 DO ji = 1, jpi 706 IF( ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0 ) THEN 707 zdz = abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) ) 708 hmld_kara(ji,jj) = gdept_n(ji,jj,ik_iso(ji,jj)) + zdz 709 ENDIF 710 END DO 711 END DO 712 713 ! If ll_found = .false. then calculate MLD using difference of zdelta_T 714 ! from the reference density/temperature 715 716 ! Prevent this section from working on land points 717 WHERE( tmask(:,:,1) /= 1.0 ) 718 ll_found = .true. 719 ENDWHERE 720 721 DO jk = 1, jpk 718 722 ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= & 719 723 & zdelta_T(:,:) 720 END DO 721 722 ! Set default value where interpolation cannot be used (ll_found=false) 723 DO jj = 1, jpj 724 DO ji = 1, jpi 724 END DO 725 726 ! Set default value where interpolation cannot be used (ll_found=false) 727 DO jj = 1, jpj 728 DO ji = 1, jpi 725 729 IF( .NOT. ll_found(ji,jj) ) & 726 & hmld_kara(ji,jj) = gdept_n(ji,jj, mbkt(ji,jj))! mbathy(ji,jj))727 END DO 728 END DO 729 730 DO jj = 1, jpj 731 DO ji = 1, jpi 732 !CDIR NOVECTOR 733 DO jk = ik_ref(ji,jj)+1, mbkt(ji,jj) !mbathy(ji,jj)734 IF( ll_found(ji,jj) ) EXIT 735 IF( ll_belowml(ji,jj,jk) ) THEN 730 & hmld_kara(ji,jj) = gdept_n(ji,jj,ikmt(ji,jj)) ! mbkt(ji,jj))! mbathy(ji,jj)) 731 END DO 732 END DO 733 734 DO jj = 1, jpj 735 DO ji = 1, jpi 736 !CDIR NOVECTOR 737 DO jk = ik_ref(ji,jj)+1, ikmt(ji,jj) !mbkt(ji,jj) !mbathy(ji,jj) 738 IF( ll_found(ji,jj) ) EXIT 739 IF( ll_belowml(ji,jj,jk) ) THEN 736 740 zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * & 737 & SIGN(1.0, zdTdz(ji,jj,jk-1) ) 738 zdT = zT_b - zT(ji,jj,jk-1) 739 zdz = zdT / zdTdz(ji,jj,jk-1) 740 hmld_kara(ji,jj) = gdept_n(ji,jj,jk-1) + zdz 741 EXIT 742 ENDIF 743 END DO 744 END DO 745 END DO 746 747 hmld_kara(:,:) = hmld_kara(:,:) * tmask(:,:,1) 748 741 & SIGN(1.0, zdTdz(ji,jj,jk-1) ) 742 zdT = zT_b - zT(ji,jj,jk-1) 743 zdz = zdT / zdTdz(ji,jj,jk-1) 744 hmld_kara(ji,jj) = gdept_n(ji,jj,jk-1) + zdz 745 EXIT 746 ENDIF 747 END DO 748 END DO 749 END DO 750 751 hmld_kara(:,:) = hmld_kara(:,:) * tmask(:,:,1) 752 749 753 IF( ln_kara_write25h ) THEN 750 754 !Double IF required as i_steps not defined if ln_kara_write25h = … … 756 760 ENDIF 757 761 ENDIF 758 759 !#if defined key_iomput 760 IF( kt /= nit000 ) THEN 761 CALL iom_put( "mldkara" , hmld_kara ) 762 763 !#if defined key_iomput 764 IF( kt /= nit000 ) THEN 765 CALL iom_put( "mldkara" , hmld_kara ) 762 766 IF( ( MOD( i_cnt_25h, 25) == 0 ) .AND. ln_kara_write25h ) & 763 767 CALL iom_put( "kara25h" , ( hmld_kara_25h / 25._wp ) ) 764 768 ENDIF 765 769 !#endif 766 770 767 771 ENDIF 768 772 ENDIF 769 770 END SUBROUTINE zdf_mxl_kara 773 774 END SUBROUTINE zdf_mxl_kara 771 775 772 776 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
