- Timestamp:
- 2017-11-17T15:40:12+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90
r8733 r8738 40 40 USE lib_fortran ! to use key_nosignedzero 41 41 #if defined key_lim3 42 USE ice , ONLY : u_ice, v_ice, jpl, pfrld, a_i_b, at_i_b 43 USE limthd_dh ! for CALL lim_thd_snwblow 44 #elif defined key_lim2 45 USE ice_2 , ONLY : u_ice, v_ice 46 USE par_ice_2 ! LIM-2 parameters 42 USE ice , ONLY : u_ice, v_ice, jpl, a_i_b, at_i_b, tm_su 43 USE icethd_dh ! for CALL ice_thd_snwblow 47 44 #endif 48 45 USE sbcblk_algo_ncar ! => turb_ncar : NCAR - CORE (Large & Yeager, 2009) … … 64 61 PUBLIC sbc_blk_init ! called in sbcmod 65 62 PUBLIC sbc_blk ! called in sbcmod 66 #if defined key_lim 2 || defined key_lim367 PUBLIC blk_ice_tau ! routine called in sbc_ice_limmodule68 PUBLIC blk_ice_flx ! routine called in sbc_ice_limmodule63 #if defined key_lim3 64 PUBLIC blk_ice_tau ! routine called in icestp module 65 PUBLIC blk_ice_flx ! routine called in icestp module 69 66 #endif 70 67 … … 96 93 REAL(wp), PARAMETER :: Ls = 2.839e6 ! latent heat of sublimation 97 94 REAL(wp), PARAMETER :: Stef = 5.67e-8 ! Stefan Boltzmann constant 98 REAL(wp), PARAMETER :: Cd_ice = 1.4e-3 ! iovi 1.63e-3 !transfer coefficient over ice95 REAL(wp), PARAMETER :: Cd_ice = 1.4e-3 ! transfer coefficient over ice 99 96 REAL(wp), PARAMETER :: albo = 0.066 ! ocean albedo assumed to be constant 100 97 ! … … 111 108 REAL(wp) :: rn_zqt ! z(q,t) : height of humidity and temperature measurements 112 109 REAL(wp) :: rn_zu ! z(u) : height of wind measurements 113 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) 114 112 ! 115 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: Cd_oce ! air-ocean drag (clem) 113 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: Cd_atm ! transfer coefficient for momentum (tau) 114 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: Ch_atm ! transfer coefficient for sensible heat (Q_sens) 115 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: Ce_atm ! tansfert coefficient for evaporation (Q_lat) 116 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_zu ! air temperature at wind speed height (needed by Lupkes 2015 bulk scheme) 117 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: q_zu ! air spec. hum. at wind speed height (needed by Lupkes 2015 bulk scheme) 118 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: cdn_oce, chn_oce, cen_oce ! needed by Lupkes 2015 bulk scheme 116 119 117 120 INTEGER :: nblk ! choice of the bulk algorithm … … 135 138 !! *** ROUTINE sbc_blk_alloc *** 136 139 !!------------------------------------------------------------------- 137 ALLOCATE( Cd_oce(jpi,jpj) , STAT=sbc_blk_alloc ) 140 ALLOCATE( Cd_atm (jpi,jpj), Ch_atm (jpi,jpj), Ce_atm (jpi,jpj), t_zu(jpi,jpj), q_zu(jpi,jpj), & 141 & cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), STAT=sbc_blk_alloc ) 138 142 ! 139 143 IF( lk_mpp ) CALL mpp_sum ( sbc_blk_alloc ) … … 167 171 & ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF, & ! bulk algorithm 168 172 & cn_dir , ln_taudif, rn_zqt, rn_zu, & 169 & rn_pfac, rn_efac, rn_vfac, ln_Cd_L12 173 & rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15 170 174 !!--------------------------------------------------------------------- 171 175 ! … … 258 262 WRITE(numout,*) ' factor applied on ocean/ice velocity rn_vfac = ', rn_vfac 259 263 WRITE(numout,*) ' (form absolute (=0) to relative winds(=1))' 264 WRITE(numout,*) ' use ice-atm drag from Lupkes2012 ln_Cd_L12 = ', ln_Cd_L12 265 WRITE(numout,*) ' use ice-atm drag from Lupkes2015 ln_Cd_L15 = ', ln_Cd_L15 260 266 ! 261 267 WRITE(numout,*) … … 364 370 REAL(wp), DIMENSION(:,:), POINTER :: zqlw, zqsb ! long wave and sensible heat fluxes 365 371 REAL(wp), DIMENSION(:,:), POINTER :: zqla, zevap ! latent heat fluxes and evaporation 366 REAL(wp), DIMENSION(:,:), POINTER :: Cd ! transfer coefficient for momentum (tau)367 REAL(wp), DIMENSION(:,:), POINTER :: Ch ! transfer coefficient for sensible heat (Q_sens)368 REAL(wp), DIMENSION(:,:), POINTER :: Ce ! tansfert coefficient for evaporation (Q_lat)369 372 REAL(wp), DIMENSION(:,:), POINTER :: zst ! surface temperature in Kelvin 370 REAL(wp), DIMENSION(:,:), POINTER :: zt_zu ! air temperature at wind speed height371 REAL(wp), DIMENSION(:,:), POINTER :: zq_zu ! air spec. hum. at wind speed height372 373 REAL(wp), DIMENSION(:,:), POINTER :: zU_zu ! bulk wind speed at height zu [m/s] 373 374 REAL(wp), DIMENSION(:,:), POINTER :: ztpot ! potential temperature of air at z=rn_zqt [K] … … 378 379 ! 379 380 CALL wrk_alloc( jpi,jpj, zwnd_i, zwnd_j, zsq, zqlw, zqsb, zqla, zevap ) 380 CALL wrk_alloc( jpi,jpj, Cd, Ch, Ce, zst, zt_zu, zq_zu ) 381 CALL wrk_alloc( jpi,jpj, zU_zu, ztpot, zrhoa ) 381 CALL wrk_alloc( jpi,jpj, zst, zU_zu, ztpot, zrhoa ) 382 382 ! 383 383 … … 426 426 zqlw(:,:) = ( sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:) ) * tmask(:,:,1) ! Long Wave 427 427 428 429 430 428 ! ----------------------------------------------------------------------------- ! 431 429 ! II Turbulent FLUXES ! … … 443 441 ! 444 442 CASE( np_NCAR ) ; CALL turb_ncar ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! NCAR-COREv2 445 & Cd, Ch, Ce, 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 ) 446 444 CASE( np_COARE_3p0 ) ; CALL turb_coare ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! COARE v3.0 447 & Cd, Ch, Ce, 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 ) 448 446 CASE( np_COARE_3p5 ) ; CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! COARE v3.5 449 & Cd, Ch, Ce, 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 ) 450 448 CASE( np_ECMWF ) ; CALL turb_ecmwf ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! ECMWF 451 & Cd, Ch, Ce, zt_zu, zq_zu, zU_zu)449 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 452 450 CASE DEFAULT 453 451 CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) … … 456 454 ! ! Compute true air density : 457 455 IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN ! At zu: (probably useless to remove zrho*grav*rn_zu from SLP...) 458 zrhoa(:,:) = rho_air( zt_zu(:,:) , zq_zu(:,:), sf(jp_slp)%fnow(:,:,1) )456 zrhoa(:,:) = rho_air( t_zu(:,:) , q_zu(:,:) , sf(jp_slp)%fnow(:,:,1) ) 459 457 ELSE ! At zt: 460 458 zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 461 459 END IF 462 460 463 Cd_oce(:,:) = Cd(:,:) ! record value of pure ocean-atm. drag (clem) 461 !! CALL iom_put( "Cd_oce", Cd_atm) ! output value of pure ocean-atm. transfer coef. 462 !! CALL iom_put( "Ch_oce", Ch_atm) ! output value of pure ocean-atm. transfer coef. 464 463 465 464 DO jj = 1, jpj ! tau module, i and j component 466 465 DO ji = 1, jpi 467 zztmp = zrhoa(ji,jj) * zU_zu(ji,jj) * Cd (ji,jj) ! using bulk wind speed466 zztmp = zrhoa(ji,jj) * zU_zu(ji,jj) * Cd_atm(ji,jj) ! using bulk wind speed 468 467 taum (ji,jj) = zztmp * wndm (ji,jj) 469 468 zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) … … 500 499 IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN 501 500 !! q_air and t_air are given at 10m (wind reference height) 502 zevap(:,:) = rn_efac*MAX( 0._wp, zqla(:,:)*Ce (:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed503 zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch (:,:)*(zst(:,:) - ztpot(:,:) ) ! Sensible Heat, using bulk wind speed501 zevap(:,:) = rn_efac*MAX( 0._wp, zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed 502 zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:) ) ! Sensible Heat, using bulk wind speed 504 503 ELSE 505 504 !! q_air and t_air are not given at 10m (wind reference height) 506 505 ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!! 507 zevap(:,:) = rn_efac*MAX( 0._wp, zqla(:,:)*Ce (:,:)*(zsq(:,:) - zq_zu(:,:) ) ) ! Evaporation !using bulk wind speed508 zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch (:,:)*(zst(:,:) - zt_zu(:,:) ) ! Sensible Heat !using bulk wind speed506 zevap(:,:) = rn_efac*MAX( 0._wp, zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - q_zu(:,:) ) ) ! Evaporation, using bulk wind speed 507 zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) ) ! Sensible Heat, using bulk wind speed 509 508 ENDIF 510 509 … … 513 512 514 513 IF(ln_ctl) THEN 515 CALL prt_ctl( tab2d_1=zqla , clinfo1=' blk_oce: zqla : ', tab2d_2=Ce , clinfo2=' Ce : ' )516 CALL prt_ctl( tab2d_1=zqsb , clinfo1=' blk_oce: zqsb : ', tab2d_2=Ch , clinfo2=' Ch: ' )514 CALL prt_ctl( tab2d_1=zqla , clinfo1=' blk_oce: zqla : ', tab2d_2=Ce_atm , clinfo2=' Ce_oce : ' ) 515 CALL prt_ctl( tab2d_1=zqsb , clinfo1=' blk_oce: zqsb : ', tab2d_2=Ch_atm , clinfo2=' Ch_oce : ' ) 517 516 CALL prt_ctl( tab2d_1=zqlw , clinfo1=' blk_oce: zqlw : ', tab2d_2=qsr, clinfo2=' qsr : ' ) 518 517 CALL prt_ctl( tab2d_1=zsq , clinfo1=' blk_oce: zsq : ', tab2d_2=zst, clinfo2=' zst : ' ) … … 566 565 ! 567 566 CALL wrk_dealloc( jpi,jpj, zwnd_i, zwnd_j, zsq, zqlw, zqsb, zqla, zevap ) 568 CALL wrk_dealloc( jpi,jpj, Cd, Ch, Ce, zst, zt_zu, zq_zu ) 569 CALL wrk_dealloc( jpi,jpj, zU_zu, ztpot, zrhoa ) 567 CALL wrk_dealloc( jpi,jpj, zst, zU_zu, ztpot, zrhoa ) 570 568 ! 571 569 IF( nn_timing == 1 ) CALL timing_stop('blk_oce') … … 573 571 END SUBROUTINE blk_oce 574 572 575 #if defined key_lim 2 || defined key_lim3573 #if defined key_lim3 576 574 577 575 SUBROUTINE blk_ice_tau … … 591 589 REAL(wp) :: zwnorm_f, zwndi_f , zwndj_f ! relative wind module and components at F-point 592 590 REAL(wp) :: zwndi_t , zwndj_t ! relative wind components at T-point 593 REAL(wp), DIMENSION(:,:), POINTER :: Cd ! transfer coefficient for momentum (tau)594 591 !!--------------------------------------------------------------------- 595 592 ! … … 597 594 ! 598 595 CALL wrk_alloc( jpi,jpj, zrhoa ) 599 CALL wrk_alloc( jpi,jpj, Cd ) 600 601 Cd(:,:) = Cd_ice 602 603 ! Make ice-atm. drag dependent on ice concentration (see Lupkes et al. 2012) (clem) 604 #if defined key_lim3 605 IF( ln_Cd_L12 ) THEN 606 CALL Cdn10_Lupkes2012( Cd ) ! calculate new drag from Lupkes(2012) equations 607 ENDIF 608 #endif 609 610 ! local scalars ( place there for vector optimisation purposes) 611 ! Computing density of air! Way denser that 1.2 over sea-ice !!! 612 !! 613 zrhoa (:,:) = rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) 614 615 !!gm brutal.... 616 utau_ice (:,:) = 0._wp 617 vtau_ice (:,:) = 0._wp 618 wndm_ice (:,:) = 0._wp 619 !!gm end 620 621 ! ----------------------------------------------------------------------------- ! 622 ! Wind components and module relative to the moving ocean ( U10m - U_ice ) ! 623 ! ----------------------------------------------------------------------------- ! 596 597 ! set transfer coefficients to default sea-ice values 598 Cd_atm(:,:) = Cd_ice 599 Ch_atm(:,:) = Cd_ice 600 Ce_atm(:,:) = Cd_ice 601 602 wndm_ice(:,:) = 0._wp !!gm brutal.... 603 604 ! ------------------------------------------------------------ ! 605 ! Wind module relative to the moving ice ( U10m - U_ice ) ! 606 ! ------------------------------------------------------------ ! 624 607 SELECT CASE( cp_ice_msh ) 625 608 CASE( 'I' ) ! B-grid ice dynamics : I-point (i.e. F-point with sea-ice indexation) … … 627 610 DO jj = 2, jpjm1 628 611 DO ji = 2, jpim1 ! B grid : NO vector opt 629 ! ... scalar wind at I-point (fld being at T-point)630 zwndi_f = 0.25 * ( sf(jp_wndi)%fnow(ji-1,jj ,1) + sf(jp_wndi)%fnow(ji ,jj ,1) &631 & + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji ,jj-1,1) ) - rn_vfac * u_ice(ji,jj)632 zwndj_f = 0.25 * ( sf(jp_wndj)%fnow(ji-1,jj ,1) + sf(jp_wndj)%fnow(ji ,jj ,1) &633 & + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji ,jj-1,1) ) - rn_vfac * v_ice(ji,jj)634 zwnorm_f = zrhoa(ji,jj) * Cd(ji,jj) * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f )635 ! ... ice stress at I-point636 utau_ice(ji,jj) = zwnorm_f * zwndi_f637 vtau_ice(ji,jj) = zwnorm_f * zwndj_f638 612 ! ... scalar wind at T-point (fld being at T-point) 639 613 zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.25 * ( u_ice(ji,jj+1) + u_ice(ji+1,jj+1) & … … 644 618 END DO 645 619 END DO 646 CALL lbc_lnk( utau_ice, 'I', -1. )647 CALL lbc_lnk( vtau_ice, 'I', -1. )648 620 CALL lbc_lnk( wndm_ice, 'T', 1. ) 649 621 ! 650 622 CASE( 'C' ) ! C-grid ice dynamics : U & V-points (same as ocean) 651 DO jj = 2, jpj 652 DO ji = fs_2, jpi! vect. opt.623 DO jj = 2, jpjm1 624 DO ji = fs_2, fs_jpim1 ! vect. opt. 653 625 zwndi_t = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) ) 654 626 zwndj_t = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji ,jj-1) + v_ice(ji,jj) ) ) … … 656 628 END DO 657 629 END DO 630 CALL lbc_lnk( wndm_ice, 'T', 1. ) 631 ! 632 END SELECT 633 634 ! Make ice-atm. drag dependent on ice concentration 635 IF ( ln_Cd_L12 ) THEN ! calculate new drag from Lupkes(2012) equations 636 CALL Cdn10_Lupkes2012( Cd_atm ) 637 Ch_atm(:,:) = Cd_atm(:,:) ! momentum and heat transfer coef. are considered identical 638 ELSEIF( ln_Cd_L15 ) THEN ! calculate new drag from Lupkes(2015) equations 639 CALL Cdn10_Lupkes2015( Cd_atm, Ch_atm ) 640 ENDIF 641 642 !! CALL iom_put( "Cd_ice", Cd_atm) ! output value of pure ice-atm. transfer coef. 643 !! CALL iom_put( "Ch_ice", Ch_atm) ! output value of pure ice-atm. transfer coef. 644 645 ! local scalars ( place there for vector optimisation purposes) 646 ! Computing density of air! Way denser that 1.2 over sea-ice !!! 647 zrhoa (:,:) = rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) 648 649 !!gm brutal.... 650 utau_ice (:,:) = 0._wp 651 vtau_ice (:,:) = 0._wp 652 !!gm end 653 654 ! ------------------------------------------------------------ ! 655 ! Wind stress relative to the moving ice ( U10m - U_ice ) ! 656 ! ------------------------------------------------------------ ! 657 SELECT CASE( cp_ice_msh ) 658 CASE( 'I' ) ! B-grid ice dynamics : I-point (i.e. F-point with sea-ice indexation) 659 DO jj = 2, jpjm1 660 DO ji = 2, jpim1 ! B grid : NO vector opt 661 ! ... scalar wind at I-point (fld being at T-point) 662 zwndi_f = 0.25 * ( sf(jp_wndi)%fnow(ji-1,jj ,1) + sf(jp_wndi)%fnow(ji ,jj ,1) & 663 & + sf(jp_wndi)%fnow(ji-1,jj-1,1) + sf(jp_wndi)%fnow(ji ,jj-1,1) ) - rn_vfac * u_ice(ji,jj) 664 zwndj_f = 0.25 * ( sf(jp_wndj)%fnow(ji-1,jj ,1) + sf(jp_wndj)%fnow(ji ,jj ,1) & 665 & + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji ,jj-1,1) ) - rn_vfac * v_ice(ji,jj) 666 ! ... ice stress at I-point 667 zwnorm_f = zrhoa(ji,jj) * Cd_atm(ji,jj) * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f ) 668 utau_ice(ji,jj) = zwnorm_f * zwndi_f 669 vtau_ice(ji,jj) = zwnorm_f * zwndj_f 670 END DO 671 END DO 672 CALL lbc_lnk( utau_ice, 'I', -1. ) 673 CALL lbc_lnk( vtau_ice, 'I', -1. ) 674 ! 675 CASE( 'C' ) ! C-grid ice dynamics : U & V-points (same as ocean) 658 676 DO jj = 2, jpjm1 659 677 DO ji = fs_2, fs_jpim1 ! vect. opt. 660 utau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd (ji,jj) * ( wndm_ice(ji+1,jj ) + wndm_ice(ji,jj) )&678 utau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji+1,jj ) + wndm_ice(ji,jj) ) & 661 679 & * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) ) 662 vtau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd (ji,jj) * ( wndm_ice(ji,jj+1 ) + wndm_ice(ji,jj) )&680 vtau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji,jj+1 ) + wndm_ice(ji,jj) ) & 663 681 & * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) ) 664 682 END DO … … 666 684 CALL lbc_lnk( utau_ice, 'U', -1. ) 667 685 CALL lbc_lnk( vtau_ice, 'V', -1. ) 668 CALL lbc_lnk( wndm_ice, 'T', 1. )669 686 ! 670 687 END SELECT … … 705 722 REAL(wp), DIMENSION(:,:) , POINTER :: zevap, zsnw ! evaporation and snw distribution after wind blowing (LIM3) 706 723 REAL(wp), DIMENSION(:,:) , POINTER :: zrhoa 707 REAL(wp), DIMENSION(:,:) , POINTER :: Cd ! transfer coefficient for momentum (tau)708 724 !!--------------------------------------------------------------------- 709 725 ! … … 711 727 ! 712 728 CALL wrk_alloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb ) 713 CALL wrk_alloc( jpi,jpj, zrhoa) 714 CALL wrk_alloc( jpi,jpj, Cd ) 715 716 Cd(:,:) = Cd_ice 717 718 ! Make ice-atm. drag dependent on ice concentration (see Lupkes et al. 2012) (clem) 719 #if defined key_lim3 720 IF( ln_Cd_L12 ) THEN 721 CALL Cdn10_Lupkes2012( Cd ) ! calculate new drag from Lupkes(2012) equations 722 ENDIF 723 #endif 724 725 ! 726 ! local scalars ( place there for vector optimisation purposes) 729 CALL wrk_alloc( jpi,jpj, zrhoa ) 730 ! 731 ! local scalars 727 732 zcoef_dqlw = 4.0 * 0.95 * Stef 728 733 zcoef_dqla = -Ls * 11637800. * (-5897.8) … … 752 757 ! ----------------------------! 753 758 754 ! ... turbulent heat fluxes 759 ! ... turbulent heat fluxes with Ch_atm recalculated in blk_ice_tau 755 760 ! Sensible Heat 756 z_qsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * C d(ji,jj) * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1))761 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)) 757 762 ! Latent Heat 758 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls * C d(ji,jj) * wndm_ice(ji,jj)&759 & * ( 11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1)) )763 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls * Ch_atm(ji,jj) * wndm_ice(ji,jj) * & 764 & ( 11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) ) 760 765 ! Latent heat sensitivity for ice (Dqla/Dt) 761 766 IF( qla_ice(ji,jj,jl) > 0._wp ) THEN 762 dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * C d(ji,jj) * wndm_ice(ji,jj) / ( zst2 ) * EXP( -5897.8 / ptsu(ji,jj,jl))767 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)) 763 768 ELSE 764 769 dqla_ice(ji,jj,jl) = 0._wp … … 766 771 767 772 ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice) 768 z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * C d(ji,jj) * wndm_ice(ji,jj)773 z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj) 769 774 770 775 ! ----------------------------! … … 786 791 CALL iom_put( 'precip' , tprecip * 86400. ) ! Total precipitation 787 792 788 #if defined key_lim3789 793 CALL wrk_alloc( jpi,jpj, zevap, zsnw ) 790 794 … … 797 801 ! --- evaporation minus precipitation --- ! 798 802 zsnw(:,:) = 0._wp 799 CALL lim_thd_snwblow( pfrld, zsnw ) ! snow distribution over ice after wind blowing800 emp_oce(:,:) = pfrld(:,:) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )803 CALL ice_thd_snwblow( (1.-at_i_b(:,:)), zsnw ) ! snow distribution over ice after wind blowing 804 emp_oce(:,:) = ( 1._wp - at_i_b(:,:) ) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw ) 801 805 emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw 802 806 emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:) 803 807 804 808 ! --- heat flux associated with emp --- ! 805 qemp_oce(:,:) = - pfrld(:,:) * zevap(:,:) * sst_m(:,:) * rcp& ! evap at sst809 qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp & ! evap at sst 806 810 & + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp & ! liquid precip at Tair 807 811 & + sprecip(:,:) * ( 1._wp - zsnw ) * & ! solid precip at min(Tair,Tsnow) … … 811 815 812 816 ! --- total solar and non solar fluxes --- ! 813 qns_tot(:,:) = pfrld(:,:) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) + qemp_ice(:,:) + qemp_oce(:,:) 814 qsr_tot(:,:) = pfrld(:,:) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 ) 817 qns_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) & 818 & + qemp_ice(:,:) + qemp_oce(:,:) 819 qsr_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 ) 815 820 816 821 ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- ! … … 824 829 825 830 CALL wrk_dealloc( jpi,jpj, zevap, zsnw ) 826 #endif827 831 828 832 !-------------------------------------------------------------------- … … 846 850 CALL wrk_dealloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb ) 847 851 CALL wrk_dealloc( jpi,jpj, zrhoa ) 848 CALL wrk_dealloc( jpi,jpj, Cd )849 852 ! 850 853 IF( nn_timing == 1 ) CALL timing_stop('blk_ice_flx') … … 973 976 974 977 #if defined key_lim3 978 975 979 SUBROUTINE Cdn10_Lupkes2012( Cd ) 976 980 !!---------------------------------------------------------------------- … … 1022 1026 1023 1027 END SUBROUTINE Cdn10_Lupkes2012 1028 1029 1030 SUBROUTINE Cdn10_Lupkes2015( Cd, Ch ) 1031 !!---------------------------------------------------------------------- 1032 !! *** ROUTINE Cdn10_Lupkes2015 *** 1033 !! 1034 !! ** pUrpose : 1lternative turbulent transfert coefficients formulation 1035 !! between sea-ice and atmosphere with distinct momentum 1036 !! and heat coefficients depending on sea-ice concentration 1037 !! and atmospheric stability (no meltponds effect for now). 1038 !! 1039 !! ** Method : The parameterization is adapted from Lupkes et al. (2015) 1040 !! and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme, 1041 !! it considers specific skin and form drags (Andreas et al. 2010) 1042 !! to compute neutral transfert coefficients for both heat and 1043 !! momemtum fluxes. Atmospheric stability effect on transfert 1044 !! coefficient is also taken into account following Louis (1979). 1045 !! 1046 !! ** References : Lupkes et al. JGR 2015 (theory) 1047 !! Lupkes et al. ECHAM6 documentation 2015 (implementation) 1048 !! 1049 !!---------------------------------------------------------------------- 1050 ! 1051 REAL(wp), DIMENSION(:,:), INTENT(inout) :: Cd 1052 REAL(wp), DIMENSION(:,:), INTENT(inout) :: Ch 1053 REAL(wp), DIMENSION(jpi,jpj) :: zst, zqo_sat, zqi_sat 1054 ! 1055 ! ECHAM6 constants 1056 REAL(wp), PARAMETER :: z0_skin_ice = 0.69e-3_wp ! Eq. 43 [m] 1057 REAL(wp), PARAMETER :: z0_form_ice = 0.57e-3_wp ! Eq. 42 [m] 1058 REAL(wp), PARAMETER :: z0_ice = 1.00e-3_wp ! Eq. 15 [m] 1059 REAL(wp), PARAMETER :: zce10 = 2.80e-3_wp ! Eq. 41 1060 REAL(wp), PARAMETER :: zbeta = 1.1_wp ! Eq. 41 1061 REAL(wp), PARAMETER :: zc = 5._wp ! Eq. 13 1062 REAL(wp), PARAMETER :: zc2 = zc * zc 1063 REAL(wp), PARAMETER :: zam = 2. * zc ! Eq. 14 1064 REAL(wp), PARAMETER :: zah = 3. * zc ! Eq. 30 1065 REAL(wp), PARAMETER :: z1_alpha = 1._wp / 0.2_wp ! Eq. 51 1066 REAL(wp), PARAMETER :: z1_alphaf = z1_alpha ! Eq. 56 1067 REAL(wp), PARAMETER :: zbetah = 1.e-3_wp ! Eq. 26 1068 REAL(wp), PARAMETER :: zgamma = 1.25_wp ! Eq. 26 1069 REAL(wp), PARAMETER :: z1_gamma = 1._wp / zgamma 1070 REAL(wp), PARAMETER :: r1_3 = 1._wp / 3._wp 1071 ! 1072 INTEGER :: ji, jj ! dummy loop indices 1073 REAL(wp) :: zthetav_os, zthetav_is, zthetav_zu 1074 REAL(wp) :: zrib_o, zrib_i 1075 REAL(wp) :: zCdn_skin_ice, zCdn_form_ice, zCdn_ice 1076 REAL(wp) :: zChn_skin_ice, zChn_form_ice 1077 REAL(wp) :: z0w, z0i, zfmi, zfmw, zfhi, zfhw 1078 REAL(wp) :: zCdn_form_tmp 1079 !!---------------------------------------------------------------------- 1080 1081 ! Momentum Neutral Transfert Coefficients (should be a constant) 1082 zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2 ! Eq. 40 1083 zCdn_skin_ice = ( vkarmn / LOG( rn_zu / z0_skin_ice + 1._wp ) )**2 ! Eq. 7 1084 zCdn_ice = zCdn_skin_ice ! Eq. 7 (cf Lupkes email for details) 1085 !zCdn_ice = 1.89e-3 ! old ECHAM5 value (cf Eq. 32) 1086 1087 ! Heat Neutral Transfert Coefficients 1088 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) 1089 1090 ! Atmospheric and Surface Variables 1091 zst(:,:) = sst_m(:,:) + rt0 ! convert SST from Celcius to Kelvin 1092 zqo_sat(:,:) = 0.98_wp * q_sat( zst(:,:) , sf(jp_slp)%fnow(:,:,1) ) ! saturation humidity over ocean [kg/kg] 1093 zqi_sat(:,:) = 0.98_wp * q_sat( tm_su(:,:), sf(jp_slp)%fnow(:,:,1) ) ! saturation humidity over ice [kg/kg] 1094 ! 1095 !! DO jj = 2, jpjm1 1096 !! DO ji = fs_2, fs_jpim1 1097 DO jj = 1, jpj 1098 DO ji = 1, jpi 1099 ! Virtual potential temperature [K] 1100 zthetav_os = zst(ji,jj) * ( 1._wp + rctv0 * zqo_sat(ji,jj) ) ! over ocean 1101 zthetav_is = tm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) ) ! ocean ice 1102 zthetav_zu = t_zu (ji,jj) * ( 1._wp + rctv0 * q_zu(ji,jj) ) ! at zu 1103 1104 ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead) 1105 zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj) )**2 ! over ocean 1106 zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2 ! over ice 1107 1108 ! Momentum and Heat Neutral Transfert Coefficients 1109 zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta ! Eq. 40 1110 zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) ) ! Eq. 53 1111 1112 ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead) 1113 z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water 1114 z0i = z0_skin_ice ! over ice (cf Lupkes email for details) 1115 IF( zrib_o <= 0._wp ) THEN 1116 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 1117 zfhw = ( 1._wp + ( zbetah * ( zthetav_os - zthetav_zu )**r1_3 / ( Chn_oce(ji,jj) * MAX(0.01, wndm(ji,jj)) ) & ! Eq. 26 1118 & )**zgamma )**z1_gamma 1119 ELSE 1120 zfmw = 1._wp / ( 1._wp + zam * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 12 1121 zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 28 1122 ENDIF 1123 1124 IF( zrib_i <= 0._wp ) THEN 1125 zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 9 1126 zfhi = 1._wp - zah * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 25 1127 ELSE 1128 zfmi = 1._wp / ( 1._wp + zam * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 11 1129 zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 27 1130 ENDIF 1131 1132 ! Momentum Transfert Coefficients (Eq. 38) 1133 Cd(ji,jj) = zCdn_skin_ice * zfmi + & 1134 & 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) ) 1135 1136 ! Heat Transfert Coefficients (Eq. 49) 1137 Ch(ji,jj) = zChn_skin_ice * zfhi + & 1138 & 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) ) 1139 ! 1140 END DO 1141 END DO 1142 !! CALL lbc_lnk_multi( Cd, 'T', 1., Ch, 'T', 1. ) 1143 ! 1144 END SUBROUTINE Cdn10_Lupkes2015 1145 1024 1146 #endif 1025 1026 1147 1027 1148 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
