- Timestamp:
- 2017-10-18T19:14:32+02:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/RK3_SRC/SBC/sbcblk.F90
r8586 r8637 40 40 USE lib_fortran ! to use key_nosignedzero 41 41 #if defined key_lim3 42 USE ice , ONLY : u_ice, v_ice, jpl, a_i_b, at_i_b 42 USE ice , ONLY : u_ice, v_ice, jpl, a_i_b, at_i_b, tm_su 43 43 USE icethd_dh ! for CALL ice_thd_snwblow 44 44 #endif … … 92 92 REAL(wp), PARAMETER :: Ls = 2.839e6 ! latent heat of sublimation 93 93 REAL(wp), PARAMETER :: Stef = 5.67e-8 ! Stefan Boltzmann constant 94 REAL(wp), PARAMETER :: Cd_ice = 1.4e-3 ! iovi 1.63e-3 !transfer coefficient over ice94 REAL(wp), PARAMETER :: Cd_ice = 1.4e-3 ! transfer coefficient over ice 95 95 REAL(wp), PARAMETER :: albo = 0.066 ! ocean albedo assumed to be constant 96 96 ! … … 108 108 REAL(wp) :: rn_zu ! z(u) : height of wind measurements 109 109 !!gm ref namelist initialize it so remove the setting to false below 110 LOGICAL :: ln_Cd_L12 = .FALSE. ! Modify the drag ice-atm and oce-atm depending on ice concentration (from Lupkes et al. JGR2012) 110 LOGICAL :: ln_Cd_L12 = .FALSE. ! Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2012) 111 LOGICAL :: ln_Cd_L15 = .FALSE. ! Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2015) 111 112 ! 112 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: Cd_oce ! air-ocean drag (clem) 113 !!cr REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: Cd_oce ! air-ocean drag (clem) 114 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: Cd_atm ! transfer coefficient for momentum (tau) 115 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: Ch_atm ! transfer coefficient for sensible heat (Q_sens) 116 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: Ce_atm ! tansfert coefficient for evaporation (Q_lat) 117 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_zu ! air temperature at wind speed height (needed by Lupkes 2015 bulk scheme) 118 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: q_zu ! air spec. hum. at wind speed height (needed by Lupkes 2015 bulk scheme) 119 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: cdn_oce, chn_oce, cen_oce ! needed by Lupkes 2015 bulk scheme 113 120 114 121 INTEGER :: nblk ! choice of the bulk algorithm … … 132 139 !! *** ROUTINE sbc_blk_alloc *** 133 140 !!------------------------------------------------------------------- 134 ALLOCATE( Cd_oce(jpi,jpj) , STAT=sbc_blk_alloc ) 141 !!cr ALLOCATE( Cd_oce(jpi,jpj) , STAT=sbc_blk_alloc ) 142 ALLOCATE( Cd_atm (jpi,jpj), Ch_atm (jpi,jpj), Ce_atm (jpi,jpj), t_zu(jpi,jpj), q_zu(jpi,jpj), & 143 & cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), STAT=sbc_blk_alloc ) 135 144 ! 136 145 IF( lk_mpp ) CALL mpp_sum ( sbc_blk_alloc ) … … 164 173 & ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF, & ! bulk algorithm 165 174 & cn_dir , ln_taudif, rn_zqt, rn_zu, & 166 & rn_pfac, rn_efac, rn_vfac, ln_Cd_L12 175 & rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15 167 176 !!--------------------------------------------------------------------- 168 177 ! … … 255 264 WRITE(numout,*) ' factor applied on ocean/ice velocity rn_vfac = ', rn_vfac 256 265 WRITE(numout,*) ' (form absolute (=0) to relative winds(=1))' 266 WRITE(numout,*) ' use ice-atm drag from Lupkes2012 ln_Cd_L12 = ', ln_Cd_L12 267 WRITE(numout,*) ' use ice-atm drag from Lupkes2015 ln_Cd_L15 = ', ln_Cd_L15 257 268 ! 258 269 WRITE(numout,*) … … 361 372 REAL(wp), DIMENSION(jpi,jpj) :: zqlw, zqsb ! long wave and sensible heat fluxes 362 373 REAL(wp), DIMENSION(jpi,jpj) :: zqla, zevap ! latent heat fluxes and evaporation 363 REAL(wp), DIMENSION(jpi,jpj) :: zCd ! transfer coefficient for momentum (tau)364 REAL(wp), DIMENSION(jpi,jpj) :: zCh ! transfer coefficient for sensible heat (Q_sens)365 REAL(wp), DIMENSION(jpi,jpj) :: zCe ! tansfert coefficient for evaporation (Q_lat)366 374 REAL(wp), DIMENSION(jpi,jpj) :: zst ! surface temperature in Kelvin 367 REAL(wp), DIMENSION(jpi,jpj) :: zt_zu ! air temperature at wind speed height368 REAL(wp), DIMENSION(jpi,jpj) :: zq_zu ! air spec. hum. at wind speed height369 375 REAL(wp), DIMENSION(jpi,jpj) :: zU_zu ! bulk wind speed at height zu [m/s] 370 376 REAL(wp), DIMENSION(jpi,jpj) :: ztpot ! potential temperature of air at z=rn_zqt [K] … … 418 424 zqlw(:,:) = ( sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:) ) * tmask(:,:,1) ! Long Wave 419 425 420 421 422 426 ! ----------------------------------------------------------------------------- ! 423 427 ! II Turbulent FLUXES ! … … 435 439 ! 436 440 CASE( np_NCAR ) ; CALL turb_ncar ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! NCAR-COREv2 437 & zCd, zCh, zCe, zt_zu, zq_zu, zU_zu)441 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 438 442 CASE( np_COARE_3p0 ) ; CALL turb_coare ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! COARE v3.0 439 & zCd, zCh, zCe, zt_zu, zq_zu, zU_zu)443 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 440 444 CASE( np_COARE_3p5 ) ; CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! COARE v3.5 441 & zCd, zCh, zCe, zt_zu, zq_zu, zU_zu)445 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 442 446 CASE( np_ECMWF ) ; CALL turb_ecmwf ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! ECMWF 443 & zCd, zCh, zCe, zt_zu, zq_zu, zU_zu)447 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 444 448 CASE DEFAULT 445 449 CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) … … 448 452 ! ! Compute true air density : 449 453 IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN ! At zu: (probably useless to remove zrho*grav*rn_zu from SLP...) 450 zrhoa(:,:) = rho_air( zt_zu(:,:) , zq_zu(:,:), sf(jp_slp)%fnow(:,:,1) )454 zrhoa(:,:) = rho_air( t_zu(:,:) , q_zu(:,:) , sf(jp_slp)%fnow(:,:,1) ) 451 455 ELSE ! At zt: 452 456 zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 453 457 END IF 454 458 455 Cd_oce(:,:) = zCd(:,:) ! record value of pure ocean-atm. drag (clem) 459 !! CALL iom_put( "Cd_oce", Cd_atm) ! output value of pure ocean-atm. transfer coef. 460 !! CALL iom_put( "Ch_oce", Ch_atm) ! output value of pure ocean-atm. transfer coef. 456 461 457 462 DO jj = 1, jpj ! tau module, i and j component 458 463 DO ji = 1, jpi 459 zztmp = zrhoa(ji,jj) * zU_zu(ji,jj) * zCd(ji,jj) ! using bulk wind speed464 zztmp = zrhoa(ji,jj) * zU_zu(ji,jj) * Cd_atm(ji,jj) ! using bulk wind speed 460 465 taum (ji,jj) = zztmp * wndm (ji,jj) 461 466 zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) … … 492 497 IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN 493 498 !! q_air and t_air are given at 10m (wind reference height) 494 zevap(:,:) = rn_efac*MAX( 0._wp, zqla(:,:)* zCe(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) )! Evaporation, using bulk wind speed495 zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)* zCh(:,:)*(zst(:,:) - ztpot(:,:)) ! Sensible Heat, using bulk wind speed499 zevap(:,:) = rn_efac*MAX( 0._wp, zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed 500 zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:) ) ! Sensible Heat, using bulk wind speed 496 501 ELSE 497 502 !! q_air and t_air are not given at 10m (wind reference height) 498 503 ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!! 499 zevap(:,:) = rn_efac*MAX( 0._wp, zqla(:,:)* zCe(:,:)*(zsq(:,:) - zq_zu(:,:) ) ) ! Evaporation !using bulk wind speed500 zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)* zCh(:,:)*(zst(:,:) - zt_zu(:,:) ) ! Sensible Heat !using bulk wind speed504 zevap(:,:) = rn_efac*MAX( 0._wp, zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - q_zu(:,:) ) ) ! Evaporation, using bulk wind speed 505 zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) ) ! Sensible Heat, using bulk wind speed 501 506 ENDIF 502 507 … … 505 510 506 511 IF(ln_ctl) THEN 507 CALL prt_ctl( tab2d_1=zqla , clinfo1=' blk_oce: zqla : ', tab2d_2= zCe , clinfo2=' Ce : ' )508 CALL prt_ctl( tab2d_1=zqsb , clinfo1=' blk_oce: zqsb : ', tab2d_2= zCh , clinfo2=' Ch: ' )512 CALL prt_ctl( tab2d_1=zqla , clinfo1=' blk_oce: zqla : ', tab2d_2=Ce_atm , clinfo2=' Ce_oce : ' ) 513 CALL prt_ctl( tab2d_1=zqsb , clinfo1=' blk_oce: zqsb : ', tab2d_2=Ch_atm , clinfo2=' Ch_oce : ' ) 509 514 CALL prt_ctl( tab2d_1=zqlw , clinfo1=' blk_oce: zqlw : ', tab2d_2=qsr, clinfo2=' qsr : ' ) 510 515 CALL prt_ctl( tab2d_1=zsq , clinfo1=' blk_oce: zsq : ', tab2d_2=zst, clinfo2=' zst : ' ) … … 574 579 !!--------------------------------------------------------------------- 575 580 INTEGER :: ji, jj ! dummy loop indices 576 REAL(wp) :: zwndi_f , zwndj_f, zwnorm_f 577 REAL(wp) :: zwndi_t , zwndj_t 578 REAL(wp), DIMENSION(jpi,jpj) :: z Cd, zrhoa! transfer coefficient for momentum (tau)581 REAL(wp) :: zwndi_f , zwndj_f, zwnorm_f ! relative wind module and components at F-point 582 REAL(wp) :: zwndi_t , zwndj_t ! relative wind components at T-point 583 REAL(wp), DIMENSION(jpi,jpj) :: zrhoa ! transfer coefficient for momentum (tau) 579 584 !!--------------------------------------------------------------------- 580 585 ! 581 586 IF( ln_timing ) CALL timing_start('blk_ice_tau') 582 587 ! 583 IF( ln_Cd_L12 ) THEN 584 CALL Cdn10_Lupkes2012( zCd ) ! air-ice drag = F(ice concentration) (see Lupkes et al., 2012) 585 ELSE 586 zCd(:,:) = Cd_ice ! constant air-ice drag 587 ENDIF 588 589 ! local scalars ( place there for vector optimisation purposes) 590 ! Computing density of air! Way denser that 1.2 over sea-ice !!! 591 !! 592 zrhoa (:,:) = rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) 593 594 !!gm brutal.... 595 utau_ice (:,:) = 0._wp 596 vtau_ice (:,:) = 0._wp 597 wndm_ice (:,:) = 0._wp 598 !!gm end 599 600 ! ----------------------------------------------------------------------------- ! 601 ! Wind components and module relative to the moving ocean ( U10m - U_ice ) ! 602 ! ----------------------------------------------------------------------------- ! 588 ! set transfer coefficients to default sea-ice values 589 Cd_atm(:,:) = Cd_ice 590 Ch_atm(:,:) = Cd_ice 591 Ce_atm(:,:) = Cd_ice 592 593 wndm_ice(:,:) = 0._wp !!gm brutal.... 594 595 ! ------------------------------------------------------------ ! 596 ! Wind module relative to the moving ice ( U10m - U_ice ) ! 597 ! ------------------------------------------------------------ ! 603 598 SELECT CASE( cp_ice_msh ) 604 599 CASE( 'I' ) ! B-grid ice dynamics : I-point (i.e. F-point with sea-ice indexation) … … 606 601 DO jj = 2, jpjm1 607 602 DO ji = 2, jpim1 ! B grid : NO vector opt 608 ! ... scalar wind at I-point (fld being at T-point)609 zwndi_f = 0.25 * ( sf(jp_wndi)%fnow(ji-1,jj ,1) + sf(jp_wndi)%fnow(ji ,jj ,1) &610 & + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji ,jj-1,1) ) - rn_vfac * u_ice(ji,jj)611 zwndj_f = 0.25 * ( sf(jp_wndj)%fnow(ji-1,jj ,1) + sf(jp_wndj)%fnow(ji ,jj ,1) &612 & + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji ,jj-1,1) ) - rn_vfac * v_ice(ji,jj)613 zwnorm_f = zrhoa(ji,jj) * zCd(ji,jj) * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f )614 ! ... ice stress at I-point615 utau_ice(ji,jj) = zwnorm_f * zwndi_f616 vtau_ice(ji,jj) = zwnorm_f * zwndj_f617 603 ! ... scalar wind at T-point (fld being at T-point) 618 604 zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.25 * ( u_ice(ji,jj+1) + u_ice(ji+1,jj+1) & … … 623 609 END DO 624 610 END DO 625 CALL lbc_lnk( utau_ice, 'I', -1. )626 CALL lbc_lnk( vtau_ice, 'I', -1. )627 611 CALL lbc_lnk( wndm_ice, 'T', 1. ) 628 612 ! 629 613 CASE( 'C' ) ! C-grid ice dynamics : U & V-points (same as ocean) 630 DO jj = 2, jpj 631 DO ji = fs_2, jpi! vect. opt.614 DO jj = 2, jpjm1 615 DO ji = fs_2, fs_jpim1 ! vect. opt. 632 616 zwndi_t = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) ) 633 617 zwndj_t = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji ,jj-1) + v_ice(ji,jj) ) ) … … 635 619 END DO 636 620 END DO 621 CALL lbc_lnk( wndm_ice, 'T', 1. ) 622 ! 623 END SELECT 624 625 ! Make ice-atm. drag dependent on ice concentration 626 IF ( ln_Cd_L12 ) THEN ! calculate new drag from Lupkes(2012) equations 627 CALL Cdn10_Lupkes2012( Cd_atm ) 628 Ch_atm(:,:) = Cd_atm(:,:) ! momentum and heat transfer coef. are considered identical 629 ELSEIF( ln_Cd_L15 ) THEN ! calculate new drag from Lupkes(2015) equations 630 CALL Cdn10_Lupkes2015( Cd_atm, Ch_atm ) 631 ENDIF 632 633 !! CALL iom_put( "Cd_ice", Cd_atm) ! output value of pure ice-atm. transfer coef. 634 !! CALL iom_put( "Ch_ice", Ch_atm) ! output value of pure ice-atm. transfer coef. 635 636 ! local scalars ( place there for vector optimisation purposes) 637 ! Computing density of air! Way denser that 1.2 over sea-ice !!! 638 zrhoa (:,:) = rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) 639 640 !!gm brutal.... 641 utau_ice (:,:) = 0._wp 642 vtau_ice (:,:) = 0._wp 643 !!gm end 644 645 ! ------------------------------------------------------------ ! 646 ! Wind stress relative to the moving ice ( U10m - U_ice ) ! 647 ! ------------------------------------------------------------ ! 648 SELECT CASE( cp_ice_msh ) 649 CASE( 'I' ) ! B-grid ice dynamics : I-point (i.e. F-point with sea-ice indexation) 650 DO jj = 2, jpjm1 651 DO ji = 2, jpim1 ! B grid : NO vector opt 652 ! ... scalar wind at I-point (fld being at T-point) 653 zwndi_f = 0.25 * ( sf(jp_wndi)%fnow(ji-1,jj ,1) + sf(jp_wndi)%fnow(ji ,jj ,1) & 654 & + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji ,jj-1,1) ) - rn_vfac * u_ice(ji,jj) 655 zwndj_f = 0.25 * ( sf(jp_wndj)%fnow(ji-1,jj ,1) + sf(jp_wndj)%fnow(ji ,jj ,1) & 656 & + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji ,jj-1,1) ) - rn_vfac * v_ice(ji,jj) 657 ! ... ice stress at I-point 658 zwnorm_f = zrhoa(ji,jj) * Cd_atm(ji,jj) * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f ) 659 utau_ice(ji,jj) = zwnorm_f * zwndi_f 660 vtau_ice(ji,jj) = zwnorm_f * zwndj_f 661 END DO 662 END DO 663 CALL lbc_lnk( utau_ice, 'I', -1. ) 664 CALL lbc_lnk( vtau_ice, 'I', -1. ) 665 ! 666 CASE( 'C' ) ! C-grid ice dynamics : U & V-points (same as ocean) 637 667 DO jj = 2, jpjm1 638 668 DO ji = fs_2, fs_jpim1 ! vect. opt. 639 utau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * zCd(ji,jj) * ( wndm_ice(ji+1,jj ) + wndm_ice(ji,jj) )&669 utau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji+1,jj ) + wndm_ice(ji,jj) ) & 640 670 & * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) ) 641 vtau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * zCd(ji,jj) * ( wndm_ice(ji,jj+1 ) + wndm_ice(ji,jj) )&671 vtau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji,jj+1 ) + wndm_ice(ji,jj) ) & 642 672 & * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) ) 643 673 END DO … … 645 675 CALL lbc_lnk( utau_ice, 'U', -1. ) 646 676 CALL lbc_lnk( vtau_ice, 'V', -1. ) 647 CALL lbc_lnk( wndm_ice, 'T', 1. )648 677 ! 649 678 END SELECT … … 684 713 REAL(wp), DIMENSION(jpi,jpj) :: zevap, zsnw ! evaporation and snw distribution after wind blowing (LIM3) 685 714 REAL(wp), DIMENSION(jpi,jpj) :: zrhoa 686 REAL(wp), DIMENSION(jpi,jpj) :: zCd ! transfer coefficient for momentum (tau)687 715 !!--------------------------------------------------------------------- 688 716 ! 689 717 IF( ln_timing ) CALL timing_start('blk_ice_flx') 690 718 ! 691 IF( ln_Cd_L12 ) THEN 692 CALL Cdn10_Lupkes2012( zCd ) ! air-ice drag = F(ice concentration) (see Lupkes et al., 2012) 693 ELSE 694 zCd(:,:) = Cd_ice ! constant air-ice drag 695 ENDIF 696 697 ! 698 ! local scalars ( place there for vector optimisation purposes) 719 ! 720 ! local scalars 699 721 zcoef_dqlw = 4.0 * 0.95 * Stef 700 722 zcoef_dqla = -Ls * 11637800. * (-5897.8) … … 724 746 ! ----------------------------! 725 747 726 ! ... turbulent heat fluxes 748 ! ... turbulent heat fluxes with Ch_atm recalculated in blk_ice_tau 727 749 ! Sensible Heat 728 z_qsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * zCd(ji,jj) * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1))750 z_qsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1)) 729 751 ! Latent Heat 730 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls * zCd(ji,jj) * wndm_ice(ji,jj)&731 & * ( 11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1)) )752 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls * Ch_atm(ji,jj) * wndm_ice(ji,jj) * & 753 & ( 11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) ) 732 754 ! Latent heat sensitivity for ice (Dqla/Dt) 733 755 IF( qla_ice(ji,jj,jl) > 0._wp ) THEN 734 dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * zCd(ji,jj) * wndm_ice(ji,jj) / ( zst2 ) * EXP( -5897.8 / ptsu(ji,jj,jl))756 dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ch_atm(ji,jj) * wndm_ice(ji,jj) / zst2 * EXP(-5897.8 / ptsu(ji,jj,jl)) 735 757 ELSE 736 758 dqla_ice(ji,jj,jl) = 0._wp … … 738 760 739 761 ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice) 740 z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * zCd(ji,jj) * wndm_ice(ji,jj)762 z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj) 741 763 742 764 ! ----------------------------! … … 935 957 #if defined key_lim3 936 958 937 SUBROUTINE Cdn10_Lupkes2012( pCd )959 SUBROUTINE Cdn10_Lupkes2012( Cd ) 938 960 !!---------------------------------------------------------------------- 939 961 !! *** ROUTINE Cdn10_Lupkes2012 *** … … 965 987 !! 966 988 !!---------------------------------------------------------------------- 967 REAL(wp), DIMENSION(:,:), INTENT( out) :: pCd ! air-ice drag coefficient989 REAL(wp), DIMENSION(:,:), INTENT(inout) :: Cd 968 990 REAL(wp), PARAMETER :: zCe = 2.23e-03_wp 969 991 REAL(wp), PARAMETER :: znu = 1._wp … … 973 995 !!---------------------------------------------------------------------- 974 996 zcoef = znu + 1._wp / ( 10._wp * zbeta ) 975 ! 997 976 998 ! generic drag over a cell partly covered by ice 977 !! pCd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) + & ! pure ocean drag978 !! & 979 !! & 999 !!Cd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) + & ! pure ocean drag 1000 !! & Cd_ice * at_i_b(:,:) + & ! pure ice drag 1001 !! & zCe * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**zmu ! change due to sea-ice morphology 980 1002 981 1003 ! ice-atm drag 982 pCd(:,:) = Cd_ice + & ! pure ice drag983 & 984 !1004 Cd(:,:) = Cd_ice + & ! pure ice drag 1005 & zCe * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp) ! change due to sea-ice morphology 1006 985 1007 END SUBROUTINE Cdn10_Lupkes2012 1008 1009 1010 SUBROUTINE Cdn10_Lupkes2015( Cd, Ch ) 1011 !!---------------------------------------------------------------------- 1012 !! *** ROUTINE Cdn10_Lupkes2015 *** 1013 !! 1014 !! ** pUrpose : 1lternative turbulent transfert coefficients formulation 1015 !! between sea-ice and atmosphere with distinct momentum 1016 !! and heat coefficients depending on sea-ice concentration 1017 !! and atmospheric stability (no meltponds effect for now). 1018 !! 1019 !! ** Method : The parameterization is adapted from Lupkes et al. (2015) 1020 !! and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme, 1021 !! it considers specific skin and form drags (Andreas et al. 2010) 1022 !! to compute neutral transfert coefficients for both heat and 1023 !! momemtum fluxes. Atmospheric stability effect on transfert 1024 !! coefficient is also taken into account following Louis (1979). 1025 !! 1026 !! ** References : Lupkes et al. JGR 2015 (theory) 1027 !! Lupkes et al. ECHAM6 documentation 2015 (implementation) 1028 !! 1029 !!---------------------------------------------------------------------- 1030 ! 1031 REAL(wp), DIMENSION(:,:), INTENT(inout) :: Cd 1032 REAL(wp), DIMENSION(:,:), INTENT(inout) :: Ch 1033 REAL(wp), DIMENSION(jpi,jpj) :: zst, zqo_sat, zqi_sat 1034 ! 1035 ! ECHAM6 constants 1036 REAL(wp), PARAMETER :: z0_skin_ice = 0.69e-3_wp ! Eq. 43 [m] 1037 REAL(wp), PARAMETER :: z0_form_ice = 0.57e-3_wp ! Eq. 42 [m] 1038 REAL(wp), PARAMETER :: z0_ice = 1.00e-3_wp ! Eq. 15 [m] 1039 REAL(wp), PARAMETER :: zce10 = 2.80e-3_wp ! Eq. 41 1040 REAL(wp), PARAMETER :: zbeta = 1.1_wp ! Eq. 41 1041 REAL(wp), PARAMETER :: zc = 5._wp ! Eq. 13 1042 REAL(wp), PARAMETER :: zc2 = zc * zc 1043 REAL(wp), PARAMETER :: zam = 2. * zc ! Eq. 14 1044 REAL(wp), PARAMETER :: zah = 3. * zc ! Eq. 30 1045 REAL(wp), PARAMETER :: z1_alpha = 1._wp / 0.2_wp ! Eq. 51 1046 REAL(wp), PARAMETER :: z1_alphaf = z1_alpha ! Eq. 56 1047 REAL(wp), PARAMETER :: zbetah = 1.e-3_wp ! Eq. 26 1048 REAL(wp), PARAMETER :: zgamma = 1.25_wp ! Eq. 26 1049 REAL(wp), PARAMETER :: z1_gamma = 1._wp / zgamma 1050 REAL(wp), PARAMETER :: r1_3 = 1._wp / 3._wp 1051 ! 1052 INTEGER :: ji, jj ! dummy loop indices 1053 REAL(wp) :: zthetav_os, zthetav_is, zthetav_zu 1054 REAL(wp) :: zrib_o, zrib_i 1055 REAL(wp) :: zCdn_skin_ice, zCdn_form_ice, zCdn_ice 1056 REAL(wp) :: zChn_skin_ice, zChn_form_ice 1057 REAL(wp) :: z0w, z0i, zfmi, zfmw, zfhi, zfhw 1058 REAL(wp) :: zCdn_form_tmp 1059 !!---------------------------------------------------------------------- 1060 1061 ! Momentum Neutral Transfert Coefficients (should be a constant) 1062 zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2 ! Eq. 40 1063 zCdn_skin_ice = ( vkarmn / LOG( rn_zu / z0_skin_ice + 1._wp ) )**2 ! Eq. 7 1064 zCdn_ice = zCdn_skin_ice ! Eq. 7 (cf Lupkes email for details) 1065 !zCdn_ice = 1.89e-3 ! old ECHAM5 value (cf Eq. 32) 1066 1067 ! Heat Neutral Transfert Coefficients 1068 zChn_skin_ice = vkarmn**2 / ( LOG( rn_zu / z0_ice + 1._wp ) * LOG( rn_zu * z1_alpha / z0_skin_ice + 1._wp ) ) ! Eq. 50 + Eq. 52 (cf Lupkes email for details) 1069 1070 ! Atmospheric and Surface Variables 1071 zst(:,:) = sst_m(:,:) + rt0 ! convert SST from Celcius to Kelvin 1072 zqo_sat(:,:) = 0.98_wp * q_sat( zst(:,:) , sf(jp_slp)%fnow(:,:,1) ) ! saturation humidity over ocean [kg/kg] 1073 zqi_sat(:,:) = 0.98_wp * q_sat( tm_su(:,:), sf(jp_slp)%fnow(:,:,1) ) ! saturation humidity over ice [kg/kg] 1074 ! 1075 DO jj = 2, jpjm1 ! reduced loop is necessary for reproducibility 1076 DO ji = fs_2, fs_jpim1 1077 ! Virtual potential temperature [K] 1078 zthetav_os = zst(ji,jj) * ( 1._wp + rctv0 * zqo_sat(ji,jj) ) ! over ocean 1079 zthetav_is = tm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) ) ! ocean ice 1080 zthetav_zu = t_zu (ji,jj) * ( 1._wp + rctv0 * q_zu(ji,jj) ) ! at zu 1081 1082 ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead) 1083 zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj) )**2 ! over ocean 1084 zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2 ! over ice 1085 1086 ! Momentum and Heat Neutral Transfert Coefficients 1087 zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta ! Eq. 40 1088 zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) ) ! Eq. 53 1089 1090 ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead) 1091 z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water 1092 z0i = z0_skin_ice ! over ice (cf Lupkes email for details) 1093 IF( zrib_o <= 0._wp ) THEN 1094 zfmw = 1._wp - zam * zrib_o / ( 1._wp + 3._wp * zc2 * Cdn_oce(ji,jj) * SQRT( -zrib_o * ( rn_zu / z0w + 1._wp ) ) ) ! Eq. 10 1095 zfhw = ( 1._wp + ( zbetah * ( zthetav_os - zthetav_zu )**r1_3 / ( Chn_oce(ji,jj) * MAX(0.01, wndm(ji,jj)) ) & ! Eq. 26 1096 & )**zgamma )**z1_gamma 1097 ELSE 1098 zfmw = 1._wp / ( 1._wp + zam * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 12 1099 zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 28 1100 ENDIF 1101 1102 IF( zrib_i <= 0._wp ) THEN 1103 zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 9 1104 zfhi = 1._wp - zah * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 25 1105 ELSE 1106 zfmi = 1._wp / ( 1._wp + zam * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 11 1107 zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 27 1108 ENDIF 1109 1110 ! Momentum Transfert Coefficients (Eq. 38) 1111 Cd(ji,jj) = zCdn_skin_ice * zfmi + & 1112 & zCdn_form_ice * ( zfmi * at_i_b(ji,jj) + zfmw * ( 1._wp - at_i_b(ji,jj) ) ) / MAX( 1.e-06, at_i_b(ji,jj) ) 1113 1114 ! Heat Transfert Coefficients (Eq. 49) 1115 Ch(ji,jj) = zChn_skin_ice * zfhi + & 1116 & zChn_form_ice * ( zfhi * at_i_b(ji,jj) + zfhw * ( 1._wp - at_i_b(ji,jj) ) ) / MAX( 1.e-06, at_i_b(ji,jj) ) 1117 ! 1118 END DO 1119 END DO 1120 CALL lbc_lnk_multi( Cd, 'T', 1., Ch, 'T', 1. ) 1121 ! 1122 END SUBROUTINE Cdn10_Lupkes2015 1123 986 1124 #endif 987 988 1125 989 1126 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.