 Timestamp:
 20200323T10:58:00+01:00 (2 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

NEMO/branches/UKMO/r12083_mixlyr_diag/src/OCE/ZDF/zdfmxl.F90
r12479 r12585 26 26 PRIVATE 27 27 28 PUBLIC zdf_mxl_tref ! called by asminc.F9029 28 PUBLIC zdf_mxl ! called by zdfphy.F90 30 29 … … 82 81 END FUNCTION zdf_mxl_alloc 83 82 84 SUBROUTINE zdf_mxl_tref()85 !!86 !! *** ROUTINE zdf_mxl_tref ***87 !!88 !! ** Purpose : Compute the mixed layer depth with temperature criteria.89 !!90 !! ** Method : The temperaturedefined mixed layer depth is required91 !! when assimilating SST in a 2D analysis.92 !!93 !! ** Action : hmld_tref94 !!95 !96 INTEGER :: ji, jj, jk ! dummy loop indices97 REAL(wp) :: t_ref ! Reference temperature98 REAL(wp) :: temp_c = 0.2 ! temperature criterion for mixed layer depth99 !!100 !101 ! Initialise array102 IF( zdf_mxl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_mxl_tref : unable to allocate arrays' )103 104 DO jj = 1, jpj105 DO ji = 1, jpi106 hmld_tref(ji,jj)=gdept_n(ji,jj,1 )107 IF(ssmask(ji,jj) > 0.)THEN108 t_ref=tsn(ji,jj,1,jp_tem)109 DO jk=2,jpk110 IF(ssmask(ji,jj)==0.)THEN111 hmld_tref(ji,jj)=gdept_n(ji,jj,jk )112 EXIT113 ELSEIF( ABS(tsn(ji,jj,jk,jp_tem)t_ref) < temp_c)THEN114 hmld_tref(ji,jj)=gdept_n(ji,jj,jk )115 ELSE116 EXIT117 ENDIF118 ENDDO119 ENDIF120 ENDDO121 ENDDO122 123 END SUBROUTINE zdf_mxl_tref124 125 83 SUBROUTINE zdf_mxl( kt ) 126 84 !! … … 210 168 END SUBROUTINE zdf_mxl 211 169 212 SUBROUTINE zdf_mxl_zint_mld( sf ) 213 !! 214 !! *** ROUTINE zdf_mxl_zint_mld *** 215 ! 216 ! Calculate verticallyinterpolated mixed layer depth diagnostic. 217 ! 218 ! This routine can calculate the mixed layer depth diagnostic suggested by 219 ! Kara et al, 2000, JGR, 105, 16803, but is more general and can calculate 220 ! verticallyinterpolated mixedlayer depth diagnostics with other parameter 221 ! settings set in the namzdf_mldzint namelist. 170 SUBROUTINE zdf_mxl_zint_mld( sf ) 171 !! 172 !! *** ROUTINE zdf_mxl_zint_mld *** 173 ! 174 ! Calculate verticallyinterpolated mixed layer depth diagnostic. 175 ! 176 ! This routine can calculate the mixed layer depth diagnostic suggested by 177 ! Kara et al, 2000, JGR, 105, 16803, but is more general and can calculate 178 ! verticallyinterpolated mixedlayer depth diagnostics with other parameter 179 ! settings set in the namzdf_mldzint namelist. 180 ! 181 ! If mld_type=1 the mixed layer depth is calculated as the depth at which the 182 ! density has increased by an amount equivalent to a temperature difference of 183 ! 0.8C at the surface. 184 ! 185 ! For other values of mld_type the mixed layer is calculated as the depth at 186 ! which the temperature differs by 0.8C from the surface temperature. 187 ! 188 ! David Acreman, Daley Calvert 189 ! 190 !! 191 192 TYPE(MXL_ZINT), INTENT(in) :: sf 193 194 ! Diagnostic criteria 195 INTEGER :: nn_mld_type ! mixed layer type 196 REAL(wp) :: rn_zref ! depth of initial T_ref 197 REAL(wp) :: rn_dT_crit ! Critical temp diff 198 REAL(wp) :: rn_iso_frac ! Fraction of rn_dT_crit used 199 200 ! Local variables 201 REAL(wp), PARAMETER :: zepsilon = 1.e30 ! local small value 202 INTEGER, DIMENSION(jpi,jpj) :: ikmt ! number of active tracer levels 203 INTEGER, DIMENSION(jpi,jpj) :: ik_ref ! index of reference level 204 INTEGER, DIMENSION(jpi,jpj) :: ik_iso ! index of last uniform temp level 205 REAL, DIMENSION(jpi,jpj,jpk) :: zT ! Temperature or density 206 REAL, DIMENSION(jpi,jpj) :: ppzdep ! depth for use in calculating d(rho) 207 REAL, DIMENSION(jpi,jpj) :: zT_ref ! reference temperature 208 REAL :: zT_b ! base temperature 209 REAL, DIMENSION(jpi,jpj,jpk) :: zdTdz ! gradient of zT 210 REAL, DIMENSION(jpi,jpj,jpk) :: zmoddT ! Absolute temperature difference 211 REAL :: zdz ! depth difference 212 REAL :: zdT ! temperature difference 213 REAL, DIMENSION(jpi,jpj) :: zdelta_T ! difference critereon 214 REAL, DIMENSION(jpi,jpj) :: zRHO1, zRHO2 ! Densities 215 INTEGER :: ji, jj, jk ! loop counter 216 217 !! 222 218 ! 223 ! If mld_type=1 the mixed layer depth is calculated as the depth at which the 224 ! density has increased by an amount equivalent to a temperature difference of 225 ! 0.8C at the surface. 226 ! 227 ! For other values of mld_type the mixed layer is calculated as the depth at 228 ! which the temperature differs by 0.8C from the surface temperature. 229 ! 230 ! David Acreman, Daley Calvert 231 ! 232 !! 233 234 TYPE(MXL_ZINT), INTENT(in) :: sf 235 236 ! Diagnostic criteria 237 INTEGER :: nn_mld_type ! mixed layer type 238 REAL(wp) :: rn_zref ! depth of initial T_ref 239 REAL(wp) :: rn_dT_crit ! Critical temp diff 240 REAL(wp) :: rn_iso_frac ! Fraction of rn_dT_crit used 241 242 ! Local variables 243 REAL(wp), PARAMETER :: zepsilon = 1.e30 ! local small value 244 INTEGER, POINTER, DIMENSION(:,:) :: ik_ref ! index of reference level 245 INTEGER, POINTER, DIMENSION(:,:) :: ik_iso ! index of last uniform temp level 246 REAL, POINTER, DIMENSION(:,:,:) :: zT ! Temperature or density 247 REAL, POINTER, DIMENSION(:,:) :: ppzdep ! depth for use in calculating d(rho) 248 REAL, POINTER, DIMENSION(:,:) :: zT_ref ! reference temperature 249 REAL :: zT_b ! base temperature 250 REAL, POINTER, DIMENSION(:,:,:) :: zdTdz ! gradient of zT 251 REAL, POINTER, DIMENSION(:,:,:) :: zmoddT ! Absolute temperature difference 252 REAL :: zdz ! depth difference 253 REAL :: zdT ! temperature difference 254 REAL, POINTER, DIMENSION(:,:) :: zdelta_T ! difference critereon 255 REAL, POINTER, DIMENSION(:,:) :: zRHO1, zRHO2 ! Densities 256 INTEGER :: ji, jj, jk ! loop counter 257 258 !! 259 ! 260 ALLOCATE( ik_ref(jpi,jpj), ik_iso(jpi,jpj), ppzdep(jpi,jpj), zT_ref(jpi,jpj), zdelta_T(jpi,jpj), & 261 zRHO1(jpi,jpj), zRHO2(jpi,jpj), zT(jpi,jpj,jpk), zdTdz(jpi,jpj,jpk), & 262 zmoddT(jpi,jpj,jpk), STAT=ji ) 263 IF( lk_mpp ) CALL mpp_sum( 'zdfmxl', ji ) 264 IF( ji /= 0 ) CALL ctl_stop( 'STOP', 'zdf_mxl_zint_mld : unable to allocate arrays' ) 265 266 ! Unpack structure 267 nn_mld_type = sf%mld_type 268 rn_zref = sf%zref 269 rn_dT_crit = sf%dT_crit 270 rn_iso_frac = sf%iso_frac 271 272 ! Set the mixed layer depth criterion at each grid point 273 IF( nn_mld_type == 0 ) THEN 274 zdelta_T(:,:) = rn_dT_crit 219 ! Unpack structure 220 nn_mld_type = sf%mld_type 221 rn_zref = sf%zref 222 rn_dT_crit = sf%dT_crit 223 rn_iso_frac = sf%iso_frac 224 225 ! Set the mixed layer depth criterion at each grid point 226 IF( nn_mld_type == 0 ) THEN 227 zdelta_T(:,:) = rn_dT_crit 228 zT(:,:,:) = rhop(:,:,:) 229 ELSE IF( nn_mld_type == 1 ) THEN 230 ppzdep(:,:)=0.0 231 call eos ( tsn(:,:,1,:), ppzdep(:,:), zRHO1(:,:) ) 232 ! Use zT temporarily as a copy of tsn with rn_dT_crit added to SST 233 ! [assumes number of tracers less than number of vertical levels] 234 zT(:,:,1:jpts)=tsn(:,:,1,1:jpts) 235 zT(:,:,jp_tem)=zT(:,:,1)+rn_dT_crit 236 CALL eos( zT(:,:,1:jpts), ppzdep(:,:), zRHO2(:,:) ) 237 zdelta_T(:,:) = abs( zRHO1(:,:)  zRHO2(:,:) ) * rau0 238 ! RHO from eos (2d version) doesn't calculate north or east halo: 239 CALL lbc_lnk( 'zdfmxl', zdelta_T, 'T', 1. ) 275 240 zT(:,:,:) = rhop(:,:,:) 276 ELSE IF( nn_mld_type == 1 ) THEN 277 ppzdep(:,:)=0.0 278 call eos ( tsn(:,:,1,:), ppzdep(:,:), zRHO1(:,:) ) 279 ! Use zT temporarily as a copy of tsn with rn_dT_crit added to SST 280 ! [assumes number of tracers less than number of vertical levels] 281 zT(:,:,1:jpts)=tsn(:,:,1,1:jpts) 282 zT(:,:,jp_tem)=zT(:,:,1)+rn_dT_crit 283 CALL eos( zT(:,:,1:jpts), ppzdep(:,:), zRHO2(:,:) ) 284 zdelta_T(:,:) = abs( zRHO1(:,:)  zRHO2(:,:) ) * rau0 285 ! RHO from eos (2d version) doesn't calculate north or east halo: 286 CALL lbc_lnk( 'zdfmlx', zdelta_T, 'T', 1. ) 287 zT(:,:,:) = rhop(:,:,:) 288 ELSE 289 zdelta_T(:,:) = rn_dT_crit 290 zT(:,:,:) = tsn(:,:,:,jp_tem) 291 END IF 292 293 ! Calculate the gradient of zT and absolute difference for use later 294 DO jk = 1 ,jpk2 295 zdTdz(:,:,jk) = ( zT(:,:,jk+1)  zT(:,:,jk) ) / e3w_n(:,:,jk+1) 296 zmoddT(:,:,jk) = abs( zT(:,:,jk+1)  zT(:,:,jk) ) 297 END DO 298 299 ! Find density/temperature at the reference level (Kara et al use 10m). 300 ! ik_ref is the index of the box centre immediately above or at the reference level 301 ! Find rn_zref in the array of model level depths and find the ref 302 ! density/temperature by linear interpolation. 303 DO jk = jpkm1, 2, 1 304 WHERE ( gdept_n(:,:,jk) > rn_zref ) 305 ik_ref(:,:) = jk  1 306 zT_ref(:,:) = zT(:,:,jk1) + zdTdz(:,:,jk1) * ( rn_zref  gdept_n(:,:,jk1) ) 307 END WHERE 308 END DO 309 310 ! If the first grid box centre is below the reference level then use the 311 ! top model level to get zT_ref 312 WHERE ( gdept_n(:,:,1) > rn_zref ) 313 zT_ref = zT(:,:,1) 314 ik_ref = 1 315 END WHERE 316 317 ! Initialize / reset 318 ll_found(:,:) = .false. 319 320 IF ( rn_iso_frac  zepsilon > 0. ) THEN 321 ! Search for a uniform density/temperature region where adjacent levels 322 ! differ by less than rn_iso_frac * deltaT. 323 ! ik_iso is the index of the last level in the uniform layer 324 ! ll_found indicates whether the mixed layer depth can be found by interpolation 325 ik_iso(:,:) = ik_ref(:,:) 326 DO jj = 1, nlcj 327 DO ji = 1, nlci 328 !CDIR NOVECTOR 329 DO jk = ik_ref(ji,jj), mbkt(ji,jj)1 330 IF ( zmoddT(ji,jj,jk) > ( rn_iso_frac * zdelta_T(ji,jj) ) ) THEN 331 ik_iso(ji,jj) = jk 332 ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) ) 333 EXIT 334 END IF 335 END DO 336 END DO 337 END DO 338 339 ! Use linear interpolation to find depth of mixed layer base where possible 340 hmld_zint(:,:) = rn_zref 341 DO jj = 1, jpj 342 DO ji = 1, jpi 343 IF (ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0) THEN 344 zdz = abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) ) 345 hmld_zint(ji,jj) = gdept_n(ji,jj,ik_iso(ji,jj)) + zdz 346 END IF 347 END DO 348 END DO 241 ELSE 242 zdelta_T(:,:) = rn_dT_crit 243 zT(:,:,:) = tsn(:,:,:,jp_tem) 349 244 END IF 350 351 ! If ll_found = .false. then calculate MLD using difference of zdelta_T 352 ! from the reference density/temperature 353 354 ! Prevent this section from working on land points 355 WHERE ( tmask(:,:,1) /= 1.0 ) 356 ll_found = .true. 357 END WHERE 358 359 DO jk=1, jpk 360 ll_belowml(:,:,jk) = abs( zT(:,:,jk)  zT_ref(:,:) ) >= zdelta_T(:,:) 361 END DO 362 363 ! Set default value where interpolation cannot be used (ll_found=false) 364 DO jj = 1, jpj 365 DO ji = 1, jpi 366 IF ( .not. ll_found(ji,jj) ) hmld_zint(ji,jj) = gdept_n(ji,jj,mbkt(ji,jj)) 367 END DO 368 END DO 369 370 DO jj = 1, jpj 371 DO ji = 1, jpi 372 !CDIR NOVECTOR 373 DO jk = ik_ref(ji,jj)+1, mbkt(ji,jj) 374 IF ( ll_found(ji,jj) ) EXIT 375 IF ( ll_belowml(ji,jj,jk) ) THEN 376 zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * SIGN(1.0, zdTdz(ji,jj,jk1) ) 377 zdT = zT_b  zT(ji,jj,jk1) 378 zdz = zdT / zdTdz(ji,jj,jk1) 379 hmld_zint(ji,jj) = gdept_n(ji,jj,jk1) + zdz 380 EXIT 381 END IF 382 END DO 383 END DO 384 END DO 385 386 hmld_zint(:,:) = hmld_zint(:,:)*tmask(:,:,1) 387 ! 388 END SUBROUTINE zdf_mxl_zint_mld 389 390 SUBROUTINE zdf_mxl_zint_htc( kt ) 391 !! 392 !! *** ROUTINE zdf_mxl_zint_htc *** 393 !! 394 !! ** Purpose : 395 !! 396 !! ** Method : 397 !! 398 399 INTEGER, INTENT(in) :: kt ! ocean timestep index 400 401 INTEGER :: ji, jj, jk 402 INTEGER :: ikmax 403 REAL(wp) :: zc, zcoef 404 ! 405 INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ilevel 406 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zthick_0, zthick 407 408 !! 409 410 IF( .NOT. ALLOCATED(ilevel) ) THEN 411 ALLOCATE( ilevel(jpi,jpj), zthick_0(jpi,jpj), & 412 & zthick(jpi,jpj), STAT=ji ) 413 IF( lk_mpp ) CALL mpp_sum( 'zdfmxl', ji ) 414 IF( ji /= 0 ) CALL ctl_stop( 'STOP', 'zdf_mxl_zint_htc : unable to allocate arrays' ) 415 ENDIF 416 417 ! Find last whole model T level above the MLD 418 ilevel(:,:) = 0 419 zthick_0(:,:) = 0._wp 420 421 DO jk = 1, jpkm1 422 DO jj = 1, jpj 423 DO ji = 1, jpi 424 zthick_0(ji,jj) = zthick_0(ji,jj) + e3t_n(ji,jj,jk) 425 IF( zthick_0(ji,jj) < hmld_zint(ji,jj) ) ilevel(ji,jj) = jk 245 246 ! Calculate the gradient of zT and absolute difference for use later 247 DO jk = 1 ,jpk2 248 zdTdz(:,:,jk) = ( zT(:,:,jk+1)  zT(:,:,jk) ) / e3w_n(:,:,jk+1) 249 zmoddT(:,:,jk) = abs( zT(:,:,jk+1)  zT(:,:,jk) ) 250 END DO 251 252 ! Find density/temperature at the reference level (Kara et al use 10m). 253 ! ik_ref is the index of the box centre immediately above or at the reference level 254 ! Find rn_zref in the array of model level depths and find the ref 255 ! density/temperature by linear interpolation. 256 DO jk = jpkm1, 2, 1 257 WHERE ( gdept_n(:,:,jk) > rn_zref ) 258 ik_ref(:,:) = jk  1 259 zT_ref(:,:) = zT(:,:,jk1) + zdTdz(:,:,jk1) * ( rn_zref  gdept_n(:,:,jk1) ) 260 END WHERE 261 END DO 262 263 ! If the first grid box centre is below the reference level then use the 264 ! top model level to get zT_ref 265 WHERE ( gdept_n(:,:,1) > rn_zref ) 266 zT_ref = zT(:,:,1) 267 ik_ref = 1 268 END WHERE 269 270 ! The number of active tracer levels is 1 less than the number of active w levels 271 ikmt(:,:) = mbkt(:,:)  1 272 273 ! Initialize / reset 274 ll_found(:,:) = .false. 275 276 IF ( rn_iso_frac  zepsilon > 0. ) THEN 277 ! Search for a uniform density/temperature region where adjacent levels 278 ! differ by less than rn_iso_frac * deltaT. 279 ! ik_iso is the index of the last level in the uniform layer 280 ! ll_found indicates whether the mixed layer depth can be found by interpolation 281 ik_iso(:,:) = ik_ref(:,:) 282 DO jj = 1, nlcj 283 DO ji = 1, nlci 284 !CDIR NOVECTOR 285 DO jk = ik_ref(ji,jj), ikmt(ji,jj)1 286 IF ( zmoddT(ji,jj,jk) > ( rn_iso_frac * zdelta_T(ji,jj) ) ) THEN 287 ik_iso(ji,jj) = jk 288 ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) ) 289 EXIT 290 END IF 291 END DO 426 292 END DO 427 293 END DO 428 WRITE(numout,*) 'zthick_0(jk =',jk,') =',zthick_0(2,2) 429 WRITE(numout,*) 'gdepw_n(jk+1 =',jk+1,') =',gdepw_n(2,2,jk+1) 294 295 ! Use linear interpolation to find depth of mixed layer base where possible 296 hmld_zint(:,:) = rn_zref 297 DO jj = 1, jpj 298 DO ji = 1, jpi 299 IF (ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0) THEN 300 zdz = abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) ) 301 hmld_zint(ji,jj) = gdept_n(ji,jj,ik_iso(ji,jj)) + zdz 302 END IF 303 END DO 304 END DO 305 END IF 306 307 ! If ll_found = .false. then calculate MLD using difference of zdelta_T 308 ! from the reference density/temperature 309 310 ! Prevent this section from working on land points 311 WHERE ( tmask(:,:,1) /= 1.0 ) 312 ll_found = .true. 313 END WHERE 314 315 DO jk=1, jpk 316 ll_belowml(:,:,jk) = abs( zT(:,:,jk)  zT_ref(:,:) ) >= zdelta_T(:,:) 430 317 END DO 431 318 432 ! Surface boundary condition 433 IF(.NOT.ln_linssh) THEN ; zthick(:,:) = 0._wp ; htc_mld(:,:) = 0._wp 434 ELSE ; zthick(:,:) = sshn(:,:) ; htc_mld(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1) 435 ENDIF 436 437 ! Deepest whole T level above the MLD 438 ikmax = MIN( MAXVAL( ilevel(:,:) ), jpkm1 ) 439 440 ! Integration down to last whole model T level 441 DO jk = 1, ikmax 442 DO jj = 1, jpj 443 DO ji = 1, jpi 444 zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, ilevel(ji,jj)  jk + 1 ) , 1 ) ) ! 0 below ilevel 445 zthick(ji,jj) = zthick(ji,jj) + zc 446 htc_mld(ji,jj) = htc_mld(ji,jj) + zc * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk) 319 ! Set default value where interpolation cannot be used (ll_found=false) 320 DO jj = 1, jpj 321 DO ji = 1, jpi 322 IF ( .not. ll_found(ji,jj) ) hmld_zint(ji,jj) = gdept_n(ji,jj,ikmt(ji,jj)) 323 END DO 324 END DO 325 326 DO jj = 1, jpj 327 DO ji = 1, jpi 328 !CDIR NOVECTOR 329 DO jk = ik_ref(ji,jj)+1, ikmt(ji,jj) 330 IF ( ll_found(ji,jj) ) EXIT 331 IF ( ll_belowml(ji,jj,jk) ) THEN 332 zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * SIGN(1.0, zdTdz(ji,jj,jk1) ) 333 zdT = zT_b  zT(ji,jj,jk1) 334 zdz = zdT / zdTdz(ji,jj,jk1) 335 hmld_zint(ji,jj) = gdept_n(ji,jj,jk1) + zdz 336 EXIT 337 END IF 447 338 END DO 448 339 END DO 449 340 END DO 450 451 ! Subsequent partial T level 452 zthick(:,:) = hmld_zint(:,:)  zthick(:,:) ! remaining thickness to reach MLD 453 454 DO jj = 1, jpj 455 DO ji = 1, jpi 456 htc_mld(ji,jj) = htc_mld(ji,jj) + tsn(ji,jj,ilevel(ji,jj)+1,jp_tem) & 457 & * MIN( e3t_n(ji,jj,ilevel(ji,jj)+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel(ji,jj)+1) 458 END DO 459 END DO 460 461 WRITE(numout,*) 'htc_mld(after) =',htc_mld(2,2) 462 463 ! Convert to heat content 464 zcoef = rau0 * rcp 465 htc_mld(:,:) = zcoef * htc_mld(:,:) 466 467 END SUBROUTINE zdf_mxl_zint_htc 468 469 SUBROUTINE zdf_mxl_zint( kt ) 470 !! 471 !! *** ROUTINE zdf_mxl_zint *** 472 !! 473 !! ** Purpose : 341 342 hmld_zint(:,:) = hmld_zint(:,:)*tmask(:,:,1) 343 ! 344 END SUBROUTINE zdf_mxl_zint_mld 345 346 SUBROUTINE zdf_mxl_zint_htc( kt ) 347 !! 348 !! *** ROUTINE zdf_mxl_zint_htc *** 474 349 !! 475 !! ** Method : 476 !! 477 478 INTEGER, INTENT(in) :: kt ! ocean timestep index 479 480 INTEGER :: ios 481 INTEGER :: jn 482 483 INTEGER :: nn_mld_diag = 0 ! number of diagnostics 484 485 INTEGER :: i_steps ! no of timesteps per hour 486 INTEGER :: ierror ! logical error message 487 488 REAL(wp) :: zdt ! timestep variable 489 490 CHARACTER(len=1) :: cmld 491 492 TYPE(MXL_ZINT) :: sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5 493 TYPE(MXL_ZINT), SAVE, DIMENSION(5) :: mld_diags 494 495 NAMELIST/namzdf_mldzint/ nn_mld_diag, sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5 496 497 !! 498 499 IF( kt == nit000 ) THEN 500 REWIND( numnam_ref ) ! Namelist namzdf_mldzint in reference namelist 501 READ ( numnam_ref, namzdf_mldzint, IOSTAT = ios, ERR = 901) 502 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in reference namelist' ) 503 504 REWIND( numnam_cfg ) ! Namelist namzdf_mldzint in configuration namelist 505 READ ( numnam_cfg, namzdf_mldzint, IOSTAT = ios, ERR = 902 ) 506 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in configuration namelist' ) 507 IF(lwm) WRITE ( numond, namzdf_mldzint ) 508 509 IF( nn_mld_diag > 5 ) CALL ctl_stop( 'STOP', 'zdf_mxl_ini: Specify no more than 5 MLD definitions' ) 510 511 mld_diags(1) = sn_mld1 512 mld_diags(2) = sn_mld2 513 mld_diags(3) = sn_mld3 514 mld_diags(4) = sn_mld4 515 mld_diags(5) = sn_mld5 516 517 IF( nn_mld_diag > 0 ) THEN 518 WRITE(numout,*) '=============== Verticallyinterpolated mixed layer ================' 519 WRITE(numout,*) '(Diagnostic number, nn_mld_type, rn_zref, rn_dT_crit, rn_iso_frac)' 520 DO jn = 1, nn_mld_diag 521 WRITE(numout,*) 'MLD criterion',jn,':' 522 WRITE(numout,*) ' nn_mld_type =', mld_diags(jn)%mld_type 523 WRITE(numout,*) ' rn_zref =' , mld_diags(jn)%zref 524 WRITE(numout,*) ' rn_dT_crit =' , mld_diags(jn)%dT_crit 525 WRITE(numout,*) ' rn_iso_frac =', mld_diags(jn)%iso_frac 526 END DO 527 WRITE(numout,*) '====================================================================' 528 ENDIF 529 ENDIF 530 531 IF( nn_mld_diag > 0 ) THEN 532 DO jn = 1, nn_mld_diag 533 WRITE(cmld,'(I1)') jn 534 IF( iom_use( "mldzint_"//cmld ) .OR. iom_use( "mldhtc_"//cmld ) ) THEN 535 CALL zdf_mxl_zint_mld( mld_diags(jn) ) 536 537 IF( iom_use( "mldzint_"//cmld ) ) THEN 538 CALL iom_put( "mldzint_"//cmld, hmld_zint(:,:) ) 539 ENDIF 540 541 IF( iom_use( "mldhtc_"//cmld ) ) THEN 542 CALL zdf_mxl_zint_htc( kt ) 543 CALL iom_put( "mldhtc_"//cmld , htc_mld(:,:) ) 544 ENDIF 545 546 IF( iom_use( "mldzint25h_"//cmld ) ) THEN 547 IF( .NOT. mld_25h_write ) mld_25h_write = .TRUE. 548 zdt = rdt 549 IF( MOD( 3600,INT(zdt) ) == 0 ) THEN 550 i_steps = 3600/INT(zdt) 551 ELSE 552 CALL ctl_stop('STOP', 'zdf_mxl_zint 25h: timestep must give MOD(3600,rdt) = 0 otherwise no hourly values are possible') 553 ENDIF 554 IF( ( mld_25h_init ) .OR. ( kt == nit000 ) ) THEN 555 i_cnt_25h = 1 556 IF( .NOT. ALLOCATED(hmld_zint_25h) ) THEN 557 ALLOCATE( hmld_zint_25h(jpi,jpj,nn_mld_diag), STAT=ierror ) 558 IF( ierror > 0 ) CALL ctl_stop( 'zdf_mxl_zint 25h: unable to allocate hmld_zint_25h' ) 559 ENDIF 560 hmld_zint_25h(:,:,jn) = hmld_zint(:,:) 561 ENDIF 562 IF( MOD( kt, i_steps ) == 0 .AND. kt .NE. nn_it000 ) THEN 563 hmld_zint_25h(:,:,jn) = hmld_zint_25h(:,:,jn) + hmld_zint(:,:) 564 ENDIF 565 IF( i_cnt_25h .EQ. 25 .AND. MOD( kt, i_steps*24) == 0 .AND. kt .NE. nn_it000 ) THEN 566 CALL iom_put( "mldzint25h_"//cmld , hmld_zint_25h(:,:,jn) / 25._wp ) 567 ENDIF 350 !! ** Purpose : 351 !! 352 !! ** Method : 353 !! 354 355 INTEGER, INTENT(in) :: kt ! ocean timestep index 356 357 INTEGER :: ji, jj, jk 358 INTEGER :: ikmax 359 REAL(wp) :: zc, zcoef 360 ! 361 INTEGER, ALLOCATABLE, DIMENSION(:,:) :: ilevel 362 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zthick_0, zthick 363 364 !! 365 366 IF( .NOT. ALLOCATED(ilevel) ) THEN 367 ALLOCATE( ilevel(jpi,jpj), zthick_0(jpi,jpj), & 368 & zthick(jpi,jpj), STAT=ji ) 369 IF( lk_mpp ) CALL mpp_sum( 'zdfmxl', ji ) 370 IF( ji /= 0 ) CALL ctl_stop( 'STOP', 'zdf_mxl_zint_htc : unable to allocate arrays' ) 371 ENDIF 372 373 ! Find last whole model T level above the MLD 374 ilevel(:,:) = 0 375 zthick_0(:,:) = 0._wp 376 377 DO jk = 1, jpkm1 378 DO jj = 1, jpj 379 DO ji = 1, jpi 380 zthick_0(ji,jj) = zthick_0(ji,jj) + e3t_n(ji,jj,jk) 381 IF( zthick_0(ji,jj) < hmld_zint(ji,jj) ) ilevel(ji,jj) = jk 382 END DO 383 END DO 384 WRITE(numout,*) 'zthick_0(jk =',jk,') =',zthick_0(2,2) 385 WRITE(numout,*) 'gdepw_n(jk+1 =',jk+1,') =',gdepw_n(2,2,jk+1) 386 END DO 387 388 ! Surface boundary condition 389 IF( ln_linssh ) THEN ; zthick(:,:) = sshn(:,:) ; htc_mld(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1) 390 ELSE ; zthick(:,:) = 0._wp ; htc_mld(:,:) = 0._wp 391 ENDIF 392 393 ! Deepest whole T level above the MLD 394 ikmax = MIN( MAXVAL( ilevel(:,:) ), jpkm1 ) 395 396 ! Integration down to last whole model T level 397 DO jk = 1, ikmax 398 DO jj = 1, jpj 399 DO ji = 1, jpi 400 zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, ilevel(ji,jj)  jk + 1 ) , 1 ) ) ! 0 below ilevel 401 zthick(ji,jj) = zthick(ji,jj) + zc 402 htc_mld(ji,jj) = htc_mld(ji,jj) + zc * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk) 403 END DO 404 END DO 405 END DO 406 407 ! Subsequent partial T level 408 zthick(:,:) = hmld_zint(:,:)  zthick(:,:) ! remaining thickness to reach MLD 409 410 DO jj = 1, jpj 411 DO ji = 1, jpi 412 htc_mld(ji,jj) = htc_mld(ji,jj) + tsn(ji,jj,ilevel(ji,jj)+1,jp_tem) & 413 & * MIN( e3t_n(ji,jj,ilevel(ji,jj)+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel(ji,jj)+1) 414 END DO 415 END DO 416 417 WRITE(numout,*) 'htc_mld(after) =',htc_mld(2,2) 418 419 ! Convert to heat content 420 zcoef = rau0 * rcp 421 htc_mld(:,:) = zcoef * htc_mld(:,:) 422 423 END SUBROUTINE zdf_mxl_zint_htc 424 425 SUBROUTINE zdf_mxl_zint( kt ) 426 !! 427 !! *** ROUTINE zdf_mxl_zint *** 428 !! 429 !! ** Purpose : 430 !! 431 !! ** Method : 432 !! 433 434 INTEGER, INTENT(in) :: kt ! ocean timestep index 435 436 INTEGER :: ios 437 INTEGER :: jn 438 439 INTEGER :: nn_mld_diag = 0 ! number of diagnostics 440 441 CHARACTER(len=1) :: cmld 442 443 TYPE(MXL_ZINT) :: sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5 444 TYPE(MXL_ZINT), SAVE, DIMENSION(5) :: mld_diags 445 446 NAMELIST/namzdf_mldzint/ nn_mld_diag, sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5 447 448 !! 449 450 IF( kt == nit000 ) THEN 451 REWIND( numnam_ref ) ! Namelist namzdf_mldzint in reference namelist 452 READ ( numnam_ref, namzdf_mldzint, IOSTAT = ios, ERR = 901) 453 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in reference namelist' ) 454 455 REWIND( numnam_cfg ) ! Namelist namzdf_mldzint in configuration namelist 456 READ ( numnam_cfg, namzdf_mldzint, IOSTAT = ios, ERR = 902 ) 457 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in configuration namelist' ) 458 IF(lwm) WRITE ( numond, namzdf_mldzint ) 459 460 IF( nn_mld_diag > 5 ) CALL ctl_stop( 'STOP', 'zdf_mxl_ini: Specify no more than 5 MLD definitions' ) 461 462 mld_diags(1) = sn_mld1 463 mld_diags(2) = sn_mld2 464 mld_diags(3) = sn_mld3 465 mld_diags(4) = sn_mld4 466 mld_diags(5) = sn_mld5 467 468 IF( lwp .AND. (nn_mld_diag > 0) ) THEN 469 WRITE(numout,*) '=============== Verticallyinterpolated mixed layer ================' 470 WRITE(numout,*) '(Diagnostic number, nn_mld_type, rn_zref, rn_dT_crit, rn_iso_frac)' 471 DO jn = 1, nn_mld_diag 472 WRITE(numout,*) 'MLD criterion',jn,':' 473 WRITE(numout,*) ' nn_mld_type =', mld_diags(jn)%mld_type 474 WRITE(numout,*) ' rn_zref =' , mld_diags(jn)%zref 475 WRITE(numout,*) ' rn_dT_crit =' , mld_diags(jn)%dT_crit 476 WRITE(numout,*) ' rn_iso_frac =', mld_diags(jn)%iso_frac 477 END DO 478 WRITE(numout,*) '====================================================================' 479 ENDIF 480 ENDIF 481 482 IF( nn_mld_diag > 0 ) THEN 483 DO jn = 1, nn_mld_diag 484 WRITE(cmld,'(I1)') jn 485 IF( iom_use( "mldzint_"//cmld ) .OR. iom_use( "mldhtc_"//cmld ) ) THEN 486 CALL zdf_mxl_zint_mld( mld_diags(jn) ) 487 488 IF( iom_use( "mldzint_"//cmld ) ) THEN 489 CALL iom_put( "mldzint_"//cmld, hmld_zint(:,:) ) 568 490 ENDIF 569 491 570 ENDIF 571 END DO 572 573 IF( mld_25h_write ) THEN 574 IF( ( MOD( kt, i_steps ) == 0 ) .OR. mld_25h_init ) THEN 575 IF (lwp) THEN 576 WRITE(numout,*) 'zdf_mxl_zint (25h) : Summed the following number of hourly values so far',i_cnt_25h 577 ENDIF 578 i_cnt_25h = i_cnt_25h + 1 579 IF( mld_25h_init ) mld_25h_init = .FALSE. 580 ENDIF 581 IF( i_cnt_25h .EQ. 25 .AND. MOD( kt, i_steps*24) == 0 .AND. kt .NE. nn_it000 ) THEN 582 i_cnt_25h = 1 583 DO jn = 1, nn_mld_diag 584 hmld_zint_25h(:,:,jn) = hmld_zint(:,:) 585 ENDDO 586 ENDIF 587 ENDIF 588 589 ENDIF 590 492 IF( iom_use( "mldhtc_"//cmld ) ) THEN 493 CALL zdf_mxl_zint_htc( kt ) 494 CALL iom_put( "mldhtc_"//cmld , htc_mld(:,:) ) 495 ENDIF 496 ENDIF 497 END DO 498 ENDIF 499 591 500 END SUBROUTINE zdf_mxl_zint 592 501
Note: See TracChangeset
for help on using the changeset viewer.