- Timestamp:
- 2015-07-08T17:13:59+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5107_hadgem3_mct/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
r5279 r5568 22 22 !! blk_oce_core : computes momentum, heat and freshwater fluxes over ocean 23 23 !! blk_ice_core : computes momentum, heat and freshwater fluxes over ice 24 !! blk_bio_meanqsr : compute daily mean short wave radiation over the ocean25 !! blk_ice_meanqsr : compute daily mean short wave radiation over the ice26 24 !! turb_core_2z : Computes turbulent transfert coefficients 27 25 !! cd_neutral_10m : Estimate of the neutral drag coefficient at 10m … … 46 44 USE sbc_ice ! Surface boundary condition: ice fields 47 45 USE lib_fortran ! to use key_nosignedzero 46 #if defined key_lim3 47 USE ice, ONLY : u_ice, v_ice, jpl, pfrld, a_i_b 48 USE limthd_dh ! for CALL lim_thd_snwblow 49 #elif defined key_lim2 50 USE ice_2, ONLY : u_ice, v_ice 51 USE par_ice_2 52 #endif 48 53 49 54 IMPLICIT NONE … … 51 56 52 57 PUBLIC sbc_blk_core ! routine called in sbcmod module 53 PUBLIC blk_ice_core ! routine called in sbc_ice_lim module 54 PUBLIC blk_ice_meanqsr ! routine called in sbc_ice_lim module 58 #if defined key_lim2 || defined key_lim3 59 PUBLIC blk_ice_core_tau ! routine called in sbc_ice_lim module 60 PUBLIC blk_ice_core_flx ! routine called in sbc_ice_lim module 61 #endif 55 62 PUBLIC turb_core_2z ! routine calles in sbcblk_mfs module 56 63 … … 195 202 ! ! compute the surface ocean fluxes using CORE bulk formulea 196 203 IF( MOD( kt - 1, nn_fsbc ) == 0 ) CALL blk_oce_core( kt, sf, sst_m, ssu_m, ssv_m ) 197 198 ! If diurnal cycle is activated, compute a daily mean short waves flux for biogeochemistery199 IF( ltrcdm2dc ) CALL blk_bio_meanqsr200 204 201 205 #if defined key_cice … … 302 306 ELSE ; qsr(:,:) = zztmp * sf(jp_qsr)%fnow(:,:,1) * tmask(:,:,1) 303 307 ENDIF 308 304 309 zqlw(:,:) = ( sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:) ) * tmask(:,:,1) ! Long Wave 305 310 ! ----------------------------------------------------------------------------- ! … … 376 381 emp (:,:) = ( zevap(:,:) & ! mass flux (evap. - precip.) 377 382 & - sf(jp_prec)%fnow(:,:,1) * rn_pfac ) * tmask(:,:,1) 378 qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:) & ! Downward Non Solar flux 383 ! 384 qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:) & ! Downward Non Solar 379 385 & - sf(jp_snow)%fnow(:,:,1) * rn_pfac * lfus & ! remove latent melting heat for solid precip 380 386 & - zevap(:,:) * pst(:,:) * rcp & ! remove evap heat content at SST … … 384 390 & * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) 385 391 ! 386 CALL iom_put( "qlw_oce", zqlw ) ! output downward longwave heat over the ocean 387 CALL iom_put( "qsb_oce", - zqsb ) ! output downward sensible heat over the ocean 388 CALL iom_put( "qla_oce", - zqla ) ! output downward latent heat over the ocean 389 CALL iom_put( "qhc_oce", qns-zqlw+zqsb+zqla ) ! output downward heat content of E-P over the ocean 390 CALL iom_put( "qns_oce", qns ) ! output downward non solar heat over the ocean 392 #if defined key_lim3 393 qns_oce(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:) ! non solar without emp (only needed by LIM3) 394 qsr_oce(:,:) = qsr(:,:) 395 #endif 396 ! 397 IF ( nn_ice == 0 ) THEN 398 CALL iom_put( "qlw_oce" , zqlw ) ! output downward longwave heat over the ocean 399 CALL iom_put( "qsb_oce" , - zqsb ) ! output downward sensible heat over the ocean 400 CALL iom_put( "qla_oce" , - zqla ) ! output downward latent heat over the ocean 401 CALL iom_put( "qemp_oce", qns-zqlw+zqsb+zqla ) ! output downward heat content of E-P over the ocean 402 CALL iom_put( "qns_oce" , qns ) ! output downward non solar heat over the ocean 403 CALL iom_put( "qsr_oce" , qsr ) ! output downward solar heat over the ocean 404 CALL iom_put( "qt_oce" , qns+qsr ) ! output total downward heat over the ocean 405 ENDIF 391 406 ! 392 407 IF(ln_ctl) THEN … … 406 421 407 422 408 SUBROUTINE blk_ice_core( pst , pui , pvi , palb , & 409 & p_taui, p_tauj, p_qns , p_qsr, & 410 & p_qla , p_dqns, p_dqla, & 411 & p_tpr , p_spr , & 412 & p_fr1 , p_fr2 , cd_grid, pdim ) 413 !!--------------------------------------------------------------------- 414 !! *** ROUTINE blk_ice_core *** 423 #if defined key_lim2 || defined key_lim3 424 SUBROUTINE blk_ice_core_tau 425 !!--------------------------------------------------------------------- 426 !! *** ROUTINE blk_ice_core_tau *** 415 427 !! 416 428 !! ** Purpose : provide the surface boundary condition over sea-ice 417 429 !! 418 !! ** Method : compute momentum, heat and freshwater exchanged 419 !! between atmosphere and sea-ice using CORE bulk 420 !! formulea, ice variables and read atmmospheric fields. 430 !! ** Method : compute momentum using CORE bulk 431 !! formulea, ice variables and read atmospheric fields. 421 432 !! NB: ice drag coefficient is assumed to be a constant 422 !! 423 !! caution : the net upward water flux has with mm/day unit 424 !!--------------------------------------------------------------------- 425 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pst ! ice surface temperature (>0, =rt0 over land) [Kelvin] 426 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pui ! ice surface velocity (i- and i- components [m/s] 427 REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pvi ! at I-point (B-grid) or U & V-point (C-grid) 428 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: palb ! ice albedo (all skies) [%] 429 REAL(wp), DIMENSION(:,:) , INTENT( out) :: p_taui ! i- & j-components of surface ice stress [N/m2] 430 REAL(wp), DIMENSION(:,:) , INTENT( out) :: p_tauj ! at I-point (B-grid) or U & V-point (C-grid) 431 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: p_qns ! non solar heat flux over ice (T-point) [W/m2] 432 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: p_qsr ! solar heat flux over ice (T-point) [W/m2] 433 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: p_qla ! latent heat flux over ice (T-point) [W/m2] 434 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: p_dqns ! non solar heat sensistivity (T-point) [W/m2] 435 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: p_dqla ! latent heat sensistivity (T-point) [W/m2] 436 REAL(wp), DIMENSION(:,:) , INTENT( out) :: p_tpr ! total precipitation (T-point) [Kg/m2/s] 437 REAL(wp), DIMENSION(:,:) , INTENT( out) :: p_spr ! solid precipitation (T-point) [Kg/m2/s] 438 REAL(wp), DIMENSION(:,:) , INTENT( out) :: p_fr1 ! 1sr fraction of qsr penetration in ice (T-point) [%] 439 REAL(wp), DIMENSION(:,:) , INTENT( out) :: p_fr2 ! 2nd fraction of qsr penetration in ice (T-point) [%] 440 CHARACTER(len=1) , INTENT(in ) :: cd_grid ! ice grid ( C or B-grid) 441 INTEGER , INTENT(in ) :: pdim ! number of ice categories 442 !! 443 INTEGER :: ji, jj, jl ! dummy loop indices 444 INTEGER :: ijpl ! number of ice categories (size of 3rd dim of input arrays) 445 REAL(wp) :: zst2, zst3 446 REAL(wp) :: zcoef_wnorm, zcoef_wnorm2, zcoef_dqlw, zcoef_dqla, zcoef_dqsb 447 REAL(wp) :: zztmp ! temporary variable 448 REAL(wp) :: zwnorm_f, zwndi_f , zwndj_f ! relative wind module and components at F-point 449 REAL(wp) :: zwndi_t , zwndj_t ! relative wind components at T-point 450 !! 451 REAL(wp), DIMENSION(:,:) , POINTER :: z_wnds_t ! wind speed ( = | U10m - U_ice | ) at T-point 452 REAL(wp), DIMENSION(:,:,:), POINTER :: z_qlw ! long wave heat flux over ice 453 REAL(wp), DIMENSION(:,:,:), POINTER :: z_qsb ! sensible heat flux over ice 454 REAL(wp), DIMENSION(:,:,:), POINTER :: z_dqlw ! long wave heat sensitivity over ice 455 REAL(wp), DIMENSION(:,:,:), POINTER :: z_dqsb ! sensible heat sensitivity over ice 456 !!--------------------------------------------------------------------- 457 ! 458 IF( nn_timing == 1 ) CALL timing_start('blk_ice_core') 459 ! 460 CALL wrk_alloc( jpi,jpj, z_wnds_t ) 461 CALL wrk_alloc( jpi,jpj,pdim, z_qlw, z_qsb, z_dqlw, z_dqsb ) 462 463 ijpl = pdim ! number of ice categories 464 433 !!--------------------------------------------------------------------- 434 INTEGER :: ji, jj ! dummy loop indices 435 REAL(wp) :: zcoef_wnorm, zcoef_wnorm2 436 REAL(wp) :: zwnorm_f, zwndi_f , zwndj_f ! relative wind module and components at F-point 437 REAL(wp) :: zwndi_t , zwndj_t ! relative wind components at T-point 438 !!--------------------------------------------------------------------- 439 ! 440 IF( nn_timing == 1 ) CALL timing_start('blk_ice_core_tau') 441 ! 465 442 ! local scalars ( place there for vector optimisation purposes) 466 443 zcoef_wnorm = rhoa * Cice 467 444 zcoef_wnorm2 = rhoa * Cice * 0.5 468 zcoef_dqlw = 4.0 * 0.95 * Stef469 zcoef_dqla = -Ls * Cice * 11637800. * (-5897.8)470 zcoef_dqsb = rhoa * cpa * Cice471 445 472 446 !!gm brutal.... 473 z_wnds_t(:,:) = 0.e0474 p_taui (:,:) = 0.e0475 p_tauj (:,:) = 0.e0447 utau_ice (:,:) = 0._wp 448 vtau_ice (:,:) = 0._wp 449 wndm_ice (:,:) = 0._wp 476 450 !!gm end 477 451 478 #if defined key_lim3479 tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) ! LIM3: make Tair available in sea-ice. WARNING allocated after call to ice_init480 #endif481 452 ! ----------------------------------------------------------------------------- ! 482 453 ! Wind components and module relative to the moving ocean ( U10m - U_ice ) ! 483 454 ! ----------------------------------------------------------------------------- ! 484 SELECT CASE( c d_grid)455 SELECT CASE( cp_ice_msh ) 485 456 CASE( 'I' ) ! B-grid ice dynamics : I-point (i.e. F-point with sea-ice indexation) 486 457 ! and scalar wind at T-point ( = | U10m - U_ice | ) (masked) … … 489 460 ! ... scalar wind at I-point (fld being at T-point) 490 461 zwndi_f = 0.25 * ( sf(jp_wndi)%fnow(ji-1,jj ,1) + sf(jp_wndi)%fnow(ji ,jj ,1) & 491 & + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji ,jj-1,1) ) - rn_vfac * pui(ji,jj)462 & + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji ,jj-1,1) ) - rn_vfac * u_ice(ji,jj) 492 463 zwndj_f = 0.25 * ( sf(jp_wndj)%fnow(ji-1,jj ,1) + sf(jp_wndj)%fnow(ji ,jj ,1) & 493 & + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji ,jj-1,1) ) - rn_vfac * pvi(ji,jj)464 & + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji ,jj-1,1) ) - rn_vfac * v_ice(ji,jj) 494 465 zwnorm_f = zcoef_wnorm * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f ) 495 466 ! ... ice stress at I-point 496 p_taui(ji,jj) = zwnorm_f * zwndi_f497 p_tauj(ji,jj) = zwnorm_f * zwndj_f467 utau_ice(ji,jj) = zwnorm_f * zwndi_f 468 vtau_ice(ji,jj) = zwnorm_f * zwndj_f 498 469 ! ... scalar wind at T-point (fld being at T-point) 499 zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.25 * ( pui(ji,jj+1) + pui(ji+1,jj+1) &500 & + pui(ji,jj ) + pui(ji+1,jj ) )501 zwndj_t = sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.25 * ( pvi(ji,jj+1) + pvi(ji+1,jj+1) &502 & + pvi(ji,jj ) + pvi(ji+1,jj ) )503 z_wnds_t(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)470 zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.25 * ( u_ice(ji,jj+1) + u_ice(ji+1,jj+1) & 471 & + u_ice(ji,jj ) + u_ice(ji+1,jj ) ) 472 zwndj_t = sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.25 * ( v_ice(ji,jj+1) + v_ice(ji+1,jj+1) & 473 & + v_ice(ji,jj ) + v_ice(ji+1,jj ) ) 474 wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 504 475 END DO 505 476 END DO 506 CALL lbc_lnk( p_taui, 'I', -1. )507 CALL lbc_lnk( p_tauj, 'I', -1. )508 CALL lbc_lnk( z_wnds_t, 'T', 1. )477 CALL lbc_lnk( utau_ice, 'I', -1. ) 478 CALL lbc_lnk( vtau_ice, 'I', -1. ) 479 CALL lbc_lnk( wndm_ice, 'T', 1. ) 509 480 ! 510 481 CASE( 'C' ) ! C-grid ice dynamics : U & V-points (same as ocean) 511 482 DO jj = 2, jpj 512 483 DO ji = fs_2, jpi ! vect. opt. 513 zwndi_t = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pui(ji-1,jj ) + pui(ji,jj) ) )514 zwndj_t = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pvi(ji ,jj-1) + pvi(ji,jj) ) )515 z_wnds_t(ji,jj)= SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)484 zwndi_t = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) ) 485 zwndj_t = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji ,jj-1) + v_ice(ji,jj) ) ) 486 wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 516 487 END DO 517 488 END DO 518 489 DO jj = 2, jpjm1 519 490 DO ji = fs_2, fs_jpim1 ! vect. opt. 520 p_taui(ji,jj) = zcoef_wnorm2 * ( z_wnds_t(ji+1,jj ) + z_wnds_t(ji,jj) ) &521 & * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * pui(ji,jj) )522 p_tauj(ji,jj) = zcoef_wnorm2 * ( z_wnds_t(ji,jj+1 ) + z_wnds_t(ji,jj) ) &523 & * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * pvi(ji,jj) )491 utau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji+1,jj ) + wndm_ice(ji,jj) ) & 492 & * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) ) 493 vtau_ice(ji,jj) = zcoef_wnorm2 * ( wndm_ice(ji,jj+1 ) + wndm_ice(ji,jj) ) & 494 & * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) ) 524 495 END DO 525 496 END DO 526 CALL lbc_lnk( p_taui, 'U', -1. )527 CALL lbc_lnk( p_tauj, 'V', -1. )528 CALL lbc_lnk( z_wnds_t, 'T', 1. )497 CALL lbc_lnk( utau_ice, 'U', -1. ) 498 CALL lbc_lnk( vtau_ice, 'V', -1. ) 499 CALL lbc_lnk( wndm_ice, 'T', 1. ) 529 500 ! 530 501 END SELECT 502 503 IF(ln_ctl) THEN 504 CALL prt_ctl(tab2d_1=utau_ice , clinfo1=' blk_ice_core: utau_ice : ', tab2d_2=vtau_ice , clinfo2=' vtau_ice : ') 505 CALL prt_ctl(tab2d_1=wndm_ice , clinfo1=' blk_ice_core: wndm_ice : ') 506 ENDIF 507 508 IF( nn_timing == 1 ) CALL timing_stop('blk_ice_core_tau') 509 510 END SUBROUTINE blk_ice_core_tau 511 512 513 SUBROUTINE blk_ice_core_flx( ptsu, palb ) 514 !!--------------------------------------------------------------------- 515 !! *** ROUTINE blk_ice_core_flx *** 516 !! 517 !! ** Purpose : provide the surface boundary condition over sea-ice 518 !! 519 !! ** Method : compute heat and freshwater exchanged 520 !! between atmosphere and sea-ice using CORE bulk 521 !! formulea, ice variables and read atmmospheric fields. 522 !! 523 !! caution : the net upward water flux has with mm/day unit 524 !!--------------------------------------------------------------------- 525 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: ptsu ! sea ice surface temperature 526 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: palb ! ice albedo (all skies) 527 !! 528 INTEGER :: ji, jj, jl ! dummy loop indices 529 REAL(wp) :: zst2, zst3 530 REAL(wp) :: zcoef_dqlw, zcoef_dqla, zcoef_dqsb 531 REAL(wp) :: zztmp, z1_lsub ! temporary variable 532 !! 533 REAL(wp), DIMENSION(:,:,:), POINTER :: z_qlw ! long wave heat flux over ice 534 REAL(wp), DIMENSION(:,:,:), POINTER :: z_qsb ! sensible heat flux over ice 535 REAL(wp), DIMENSION(:,:,:), POINTER :: z_dqlw ! long wave heat sensitivity over ice 536 REAL(wp), DIMENSION(:,:,:), POINTER :: z_dqsb ! sensible heat sensitivity over ice 537 REAL(wp), DIMENSION(:,:) , POINTER :: zevap, zsnw ! evaporation and snw distribution after wind blowing (LIM3) 538 !!--------------------------------------------------------------------- 539 ! 540 IF( nn_timing == 1 ) CALL timing_start('blk_ice_core_flx') 541 ! 542 CALL wrk_alloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb ) 543 544 ! local scalars ( place there for vector optimisation purposes) 545 zcoef_dqlw = 4.0 * 0.95 * Stef 546 zcoef_dqla = -Ls * Cice * 11637800. * (-5897.8) 547 zcoef_dqsb = rhoa * cpa * Cice 531 548 532 549 zztmp = 1. / ( 1. - albo ) 533 550 ! ! ========================== ! 534 DO jl = 1, ijpl! Loop over ice categories !551 DO jl = 1, jpl ! Loop over ice categories ! 535 552 ! ! ========================== ! 536 553 DO jj = 1 , jpj … … 539 556 ! I Radiative FLUXES ! 540 557 ! ----------------------------! 541 zst2 = p st(ji,jj,jl) * pst(ji,jj,jl)542 zst3 = p st(ji,jj,jl) * zst2558 zst2 = ptsu(ji,jj,jl) * ptsu(ji,jj,jl) 559 zst3 = ptsu(ji,jj,jl) * zst2 543 560 ! Short Wave (sw) 544 p_qsr(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj)561 qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj) 545 562 ! Long Wave (lw) 546 z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * p st(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)563 z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 547 564 ! lw sensitivity 548 565 z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3 … … 554 571 ! ... turbulent heat fluxes 555 572 ! Sensible Heat 556 z_qsb(ji,jj,jl) = rhoa * cpa * Cice * z_wnds_t(ji,jj) * ( pst(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) )573 z_qsb(ji,jj,jl) = rhoa * cpa * Cice * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) ) 557 574 ! Latent Heat 558 p_qla(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls * Cice * z_wnds_t(ji,jj) &559 & * ( 11637800. * EXP( -5897.8 / p st(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1) ) )560 561 IF( p_qla(ji,jj,jl) > 0._wp ) THEN562 p_dqla(ji,jj,jl) = rn_efac * zcoef_dqla * z_wnds_t(ji,jj) / ( zst2 ) * EXP( -5897.8 / pst(ji,jj,jl) )575 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa * Ls * Cice * wndm_ice(ji,jj) & 576 & * ( 11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / rhoa - sf(jp_humi)%fnow(ji,jj,1) ) ) 577 ! Latent heat sensitivity for ice (Dqla/Dt) 578 IF( qla_ice(ji,jj,jl) > 0._wp ) THEN 579 dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * wndm_ice(ji,jj) / ( zst2 ) * EXP( -5897.8 / ptsu(ji,jj,jl) ) 563 580 ELSE 564 p_dqla(ji,jj,jl) = 0._wp581 dqla_ice(ji,jj,jl) = 0._wp 565 582 ENDIF 566 583 567 584 ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice) 568 z_dqsb(ji,jj,jl) = zcoef_dqsb * z_wnds_t(ji,jj)585 z_dqsb(ji,jj,jl) = zcoef_dqsb * wndm_ice(ji,jj) 569 586 570 587 ! ----------------------------! … … 572 589 ! ----------------------------! 573 590 ! Downward Non Solar flux 574 p_qns (ji,jj,jl) = z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - p_qla(ji,jj,jl)591 qns_ice (ji,jj,jl) = z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl) 575 592 ! Total non solar heat flux sensitivity for ice 576 p_dqns(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + p_dqla(ji,jj,jl) )593 dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) ) 577 594 END DO 578 595 ! … … 581 598 END DO 582 599 ! 600 tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac ! total precipitation [kg/m2/s] 601 sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac ! solid precipitation [kg/m2/s] 602 CALL iom_put( 'snowpre', sprecip * 86400. ) ! Snow precipitation 603 CALL iom_put( 'precip' , tprecip * 86400. ) ! Total precipitation 604 605 #if defined key_lim3 606 CALL wrk_alloc( jpi,jpj, zevap, zsnw ) 607 608 ! --- evaporation --- ! 609 z1_lsub = 1._wp / Lsub 610 evap_ice (:,:,:) = qla_ice (:,:,:) * z1_lsub ! sublimation 611 devap_ice(:,:,:) = dqla_ice(:,:,:) * z1_lsub 612 zevap (:,:) = emp(:,:) + tprecip(:,:) ! evaporation over ocean 613 614 ! --- evaporation minus precipitation --- ! 615 zsnw(:,:) = 0._wp 616 CALL lim_thd_snwblow( pfrld, zsnw ) ! snow distribution over ice after wind blowing 617 emp_oce(:,:) = pfrld(:,:) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw ) 618 emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw 619 emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:) 620 621 ! --- heat flux associated with emp --- ! 622 qemp_oce(:,:) = - pfrld(:,:) * zevap(:,:) * sst_m(:,:) * rcp & ! evap at sst 623 & + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp & ! liquid precip at Tair 624 & + sprecip(:,:) * ( 1._wp - zsnw ) * & ! solid precip at min(Tair,Tsnow) 625 & ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus ) 626 qemp_ice(:,:) = sprecip(:,:) * zsnw * & ! solid precip (only) 627 & ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus ) 628 629 ! --- total solar and non solar fluxes --- ! 630 qns_tot(:,:) = pfrld(:,:) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) + qemp_ice(:,:) + qemp_oce(:,:) 631 qsr_tot(:,:) = pfrld(:,:) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 ) 632 633 ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- ! 634 qprec_ice(:,:) = rhosn * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0_snow ) - rt0 ) * cpic * tmask(:,:,1) - lfus ) 635 636 CALL wrk_dealloc( jpi,jpj, zevap, zsnw ) 637 #endif 638 583 639 !-------------------------------------------------------------------- 584 640 ! FRACTIONs of net shortwave radiation which is not absorbed in the … … 586 642 ! ( Maykut and Untersteiner, 1971 ; Ebert and Curry, 1993 ) 587 643 ! 588 p_fr1(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 589 p_fr2(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 590 ! 591 p_tpr(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac ! total precipitation [kg/m2/s] 592 p_spr(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac ! solid precipitation [kg/m2/s] 593 CALL iom_put( 'snowpre', p_spr * 86400. ) ! Snow precipitation 594 CALL iom_put( 'precip' , p_tpr * 86400. ) ! Total precipitation 644 fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 645 fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 646 ! 595 647 ! 596 648 IF(ln_ctl) THEN 597 CALL prt_ctl(tab3d_1=p_qla , clinfo1=' blk_ice_core: p_qla : ', tab3d_2=z_qsb , clinfo2=' z_qsb : ', kdim=ijpl) 598 CALL prt_ctl(tab3d_1=z_qlw , clinfo1=' blk_ice_core: z_qlw : ', tab3d_2=p_dqla , clinfo2=' p_dqla : ', kdim=ijpl) 599 CALL prt_ctl(tab3d_1=z_dqsb , clinfo1=' blk_ice_core: z_dqsb : ', tab3d_2=z_dqlw , clinfo2=' z_dqlw : ', kdim=ijpl) 600 CALL prt_ctl(tab3d_1=p_dqns , clinfo1=' blk_ice_core: p_dqns : ', tab3d_2=p_qsr , clinfo2=' p_qsr : ', kdim=ijpl) 601 CALL prt_ctl(tab3d_1=pst , clinfo1=' blk_ice_core: pst : ', tab3d_2=p_qns , clinfo2=' p_qns : ', kdim=ijpl) 602 CALL prt_ctl(tab2d_1=p_tpr , clinfo1=' blk_ice_core: p_tpr : ', tab2d_2=p_spr , clinfo2=' p_spr : ') 603 CALL prt_ctl(tab2d_1=p_taui , clinfo1=' blk_ice_core: p_taui : ', tab2d_2=p_tauj , clinfo2=' p_tauj : ') 604 CALL prt_ctl(tab2d_1=z_wnds_t, clinfo1=' blk_ice_core: z_wnds_t : ') 605 ENDIF 606 607 CALL wrk_dealloc( jpi,jpj, z_wnds_t ) 608 CALL wrk_dealloc( jpi,jpj, pdim, z_qlw, z_qsb, z_dqlw, z_dqsb ) 609 ! 610 IF( nn_timing == 1 ) CALL timing_stop('blk_ice_core') 611 ! 612 END SUBROUTINE blk_ice_core 613 614 615 SUBROUTINE blk_bio_meanqsr 616 !!--------------------------------------------------------------------- 617 !! *** ROUTINE blk_bio_meanqsr 618 !! 619 !! ** Purpose : provide daily qsr_mean for PISCES when 620 !! analytic diurnal cycle is applied in physic 621 !! 622 !! ** Method : add part where there is no ice 623 !! 624 !!--------------------------------------------------------------------- 625 IF( nn_timing == 1 ) CALL timing_start('blk_bio_meanqsr') 626 ! 627 qsr_mean(:,:) = (1. - albo ) * sf(jp_qsr)%fnow(:,:,1) 628 ! 629 IF( nn_timing == 1 ) CALL timing_stop('blk_bio_meanqsr') 630 ! 631 END SUBROUTINE blk_bio_meanqsr 632 633 634 SUBROUTINE blk_ice_meanqsr( palb, p_qsr_mean, pdim ) 635 !!--------------------------------------------------------------------- 636 !! 637 !! ** Purpose : provide the daily qsr_mean over sea_ice for PISCES when 638 !! analytic diurnal cycle is applied in physic 639 !! 640 !! ** Method : compute qsr 641 !! 642 !!--------------------------------------------------------------------- 643 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: palb ! ice albedo (clear sky) (alb_ice_cs) [%] 644 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: p_qsr_mean ! solar heat flux over ice (T-point) [W/m2] 645 INTEGER , INTENT(in ) :: pdim ! number of ice categories 646 ! 647 INTEGER :: ijpl ! number of ice categories (size of 3rd dim of input arrays) 648 INTEGER :: ji, jj, jl ! dummy loop indices 649 REAL(wp) :: zztmp ! temporary variable 650 !!--------------------------------------------------------------------- 651 IF( nn_timing == 1 ) CALL timing_start('blk_ice_meanqsr') 652 ! 653 ijpl = pdim ! number of ice categories 654 zztmp = 1. / ( 1. - albo ) 655 ! ! ========================== ! 656 DO jl = 1, ijpl ! Loop over ice categories ! 657 ! ! ========================== ! 658 DO jj = 1 , jpj 659 DO ji = 1, jpi 660 p_qsr_mean(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr_mean(ji,jj) 661 END DO 662 END DO 663 END DO 664 ! 665 IF( nn_timing == 1 ) CALL timing_stop('blk_ice_meanqsr') 666 ! 667 END SUBROUTINE blk_ice_meanqsr 668 649 CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice_core: qla_ice : ', tab3d_2=z_qsb , clinfo2=' z_qsb : ', kdim=jpl) 650 CALL prt_ctl(tab3d_1=z_qlw , clinfo1=' blk_ice_core: z_qlw : ', tab3d_2=dqla_ice, clinfo2=' dqla_ice : ', kdim=jpl) 651 CALL prt_ctl(tab3d_1=z_dqsb , clinfo1=' blk_ice_core: z_dqsb : ', tab3d_2=z_dqlw , clinfo2=' z_dqlw : ', kdim=jpl) 652 CALL prt_ctl(tab3d_1=dqns_ice, clinfo1=' blk_ice_core: dqns_ice : ', tab3d_2=qsr_ice , clinfo2=' qsr_ice : ', kdim=jpl) 653 CALL prt_ctl(tab3d_1=ptsu , clinfo1=' blk_ice_core: ptsu : ', tab3d_2=qns_ice , clinfo2=' qns_ice : ', kdim=jpl) 654 CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice_core: tprecip : ', tab2d_2=sprecip , clinfo2=' sprecip : ') 655 ENDIF 656 657 CALL wrk_dealloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb ) 658 ! 659 IF( nn_timing == 1 ) CALL timing_stop('blk_ice_core_flx') 660 661 END SUBROUTINE blk_ice_core_flx 662 #endif 669 663 670 664 SUBROUTINE turb_core_2z( zt, zu, sst, T_zt, q_sat, q_zt, dU, &
Note: See TracChangeset
for help on using the changeset viewer.