- Timestamp:
- 2017-12-01T14:53:57+01:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90
r8877 r8879 17 17 !! ==> based on AeroBulk (http://aerobulk.sourceforge.net/) 18 18 !! 4.0 ! 2016-10 (G. Madec) introduce a sbc_blk_init routine 19 !! 4.0 ! 2016-10 (M. Vancoppenolle) Introduce Jules emulator (M. Vancoppenolle) 19 20 !!---------------------------------------------------------------------- 20 21 … … 40 41 USE lib_fortran ! to use key_nosignedzero 41 42 #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 43 USE ice , ONLY : u_ice, v_ice, jpl, a_i_b, at_i_b, tm_su 44 USE icethd_dh ! for CALL ice_thd_snwblow 45 USE icethd_zdf, ONLY: rn_cnd_s ! for blk_ice_qcn 47 46 #endif 48 47 USE sbcblk_algo_ncar ! => turb_ncar : NCAR - CORE (Large & Yeager, 2009) … … 64 63 PUBLIC sbc_blk_init ! called in sbcmod 65 64 PUBLIC sbc_blk ! called in sbcmod 66 #if defined key_lim2 || defined key_lim3 67 PUBLIC blk_ice_tau ! routine called in sbc_ice_lim module 68 PUBLIC blk_ice_flx ! routine called in sbc_ice_lim module 69 #endif 65 #if defined key_lim3 66 PUBLIC blk_ice_tau ! routine called in icestp module 67 PUBLIC blk_ice_flx ! routine called in icestp module 68 PUBLIC blk_ice_qcn ! routine called in icestp module 69 #endif 70 70 71 71 !!Lolo: should ultimately be moved in the module with all physical constants ? … … 96 96 REAL(wp), PARAMETER :: Ls = 2.839e6 ! latent heat of sublimation 97 97 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 ice98 REAL(wp), PARAMETER :: Cd_ice = 1.4e-3 ! transfer coefficient over ice 99 99 REAL(wp), PARAMETER :: albo = 0.066 ! ocean albedo assumed to be constant 100 100 ! … … 111 111 REAL(wp) :: rn_zqt ! z(q,t) : height of humidity and temperature measurements 112 112 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) 113 LOGICAL :: ln_Cd_L12 = .FALSE. ! Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2012) 114 LOGICAL :: ln_Cd_L15 = .FALSE. ! Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2015) 114 115 ! 115 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: Cd_oce ! air-ocean drag (clem) 116 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: Cd_atm ! transfer coefficient for momentum (tau) 117 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: Ch_atm ! transfer coefficient for sensible heat (Q_sens) 118 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: Ce_atm ! tansfert coefficient for evaporation (Q_lat) 119 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_zu ! air temperature at wind speed height (needed by Lupkes 2015 bulk scheme) 120 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: q_zu ! air spec. hum. at wind speed height (needed by Lupkes 2015 bulk scheme) 121 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: cdn_oce, chn_oce, cen_oce ! needed by Lupkes 2015 bulk scheme 116 122 117 123 INTEGER :: nblk ! choice of the bulk algorithm … … 135 141 !! *** ROUTINE sbc_blk_alloc *** 136 142 !!------------------------------------------------------------------- 137 ALLOCATE( Cd_oce(jpi,jpj) , STAT=sbc_blk_alloc ) 143 ALLOCATE( Cd_atm (jpi,jpj), Ch_atm (jpi,jpj), Ce_atm (jpi,jpj), t_zu(jpi,jpj), q_zu(jpi,jpj), & 144 & cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), STAT=sbc_blk_alloc ) 138 145 ! 139 146 IF( lk_mpp ) CALL mpp_sum ( sbc_blk_alloc ) … … 167 174 & ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF, & ! bulk algorithm 168 175 & cn_dir , ln_taudif, rn_zqt, rn_zu, & 169 & rn_pfac, rn_efac, rn_vfac, ln_Cd_L12 176 & rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15 170 177 !!--------------------------------------------------------------------- 171 178 ! … … 258 265 WRITE(numout,*) ' factor applied on ocean/ice velocity rn_vfac = ', rn_vfac 259 266 WRITE(numout,*) ' (form absolute (=0) to relative winds(=1))' 267 WRITE(numout,*) ' use ice-atm drag from Lupkes2012 ln_Cd_L12 = ', ln_Cd_L12 268 WRITE(numout,*) ' use ice-atm drag from Lupkes2015 ln_Cd_L15 = ', ln_Cd_L15 260 269 ! 261 270 WRITE(numout,*) … … 364 373 REAL(wp), DIMENSION(:,:), POINTER :: zqlw, zqsb ! long wave and sensible heat fluxes 365 374 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 375 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 376 REAL(wp), DIMENSION(:,:), POINTER :: zU_zu ! bulk wind speed at height zu [m/s] 373 377 REAL(wp), DIMENSION(:,:), POINTER :: ztpot ! potential temperature of air at z=rn_zqt [K] … … 378 382 ! 379 383 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 ) 384 CALL wrk_alloc( jpi,jpj, zst, zU_zu, ztpot, zrhoa ) 382 385 ! 383 386 … … 426 429 zqlw(:,:) = ( sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:) ) * tmask(:,:,1) ! Long Wave 427 430 428 429 430 431 ! ----------------------------------------------------------------------------- ! 431 432 ! II Turbulent FLUXES ! … … 443 444 ! 444 445 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)446 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 446 447 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)448 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 448 449 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)450 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 450 451 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)452 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 452 453 CASE DEFAULT 453 454 CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) … … 456 457 ! ! Compute true air density : 457 458 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) )459 zrhoa(:,:) = rho_air( t_zu(:,:) , q_zu(:,:) , sf(jp_slp)%fnow(:,:,1) ) 459 460 ELSE ! At zt: 460 461 zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 461 462 END IF 462 463 463 Cd_oce(:,:) = Cd(:,:) ! record value of pure ocean-atm. drag (clem) 464 !! CALL iom_put( "Cd_oce", Cd_atm) ! output value of pure ocean-atm. transfer coef. 465 !! CALL iom_put( "Ch_oce", Ch_atm) ! output value of pure ocean-atm. transfer coef. 464 466 465 467 DO jj = 1, jpj ! tau module, i and j component 466 468 DO ji = 1, jpi 467 zztmp = zrhoa(ji,jj) * zU_zu(ji,jj) * Cd (ji,jj) ! using bulk wind speed469 zztmp = zrhoa(ji,jj) * zU_zu(ji,jj) * Cd_atm(ji,jj) ! using bulk wind speed 468 470 taum (ji,jj) = zztmp * wndm (ji,jj) 469 471 zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) … … 500 502 IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN 501 503 !! 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 speed504 zevap(:,:) = rn_efac*MAX( 0._wp, zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed 505 zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:) ) ! Sensible Heat, using bulk wind speed 504 506 ELSE 505 507 !! q_air and t_air are not given at 10m (wind reference height) 506 508 ! 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 speed509 zevap(:,:) = rn_efac*MAX( 0._wp, zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - q_zu(:,:) ) ) ! Evaporation, using bulk wind speed 510 zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) ) ! Sensible Heat, using bulk wind speed 509 511 ENDIF 510 512 … … 513 515 514 516 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: ' )517 CALL prt_ctl( tab2d_1=zqla , clinfo1=' blk_oce: zqla : ', tab2d_2=Ce_atm , clinfo2=' Ce_oce : ' ) 518 CALL prt_ctl( tab2d_1=zqsb , clinfo1=' blk_oce: zqsb : ', tab2d_2=Ch_atm , clinfo2=' Ch_oce : ' ) 517 519 CALL prt_ctl( tab2d_1=zqlw , clinfo1=' blk_oce: zqlw : ', tab2d_2=qsr, clinfo2=' qsr : ' ) 518 520 CALL prt_ctl( tab2d_1=zsq , clinfo1=' blk_oce: zsq : ', tab2d_2=zst, clinfo2=' zst : ' ) … … 566 568 ! 567 569 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 ) 570 CALL wrk_dealloc( jpi,jpj, zst, zU_zu, ztpot, zrhoa ) 570 571 ! 571 572 IF( nn_timing == 1 ) CALL timing_stop('blk_oce') … … 573 574 END SUBROUTINE blk_oce 574 575 575 #if defined key_lim 2 || defined key_lim3576 #if defined key_lim3 576 577 577 578 SUBROUTINE blk_ice_tau … … 591 592 REAL(wp) :: zwnorm_f, zwndi_f , zwndj_f ! relative wind module and components at F-point 592 593 REAL(wp) :: zwndi_t , zwndj_t ! relative wind components at T-point 593 REAL(wp), DIMENSION(:,:), POINTER :: Cd ! transfer coefficient for momentum (tau)594 594 !!--------------------------------------------------------------------- 595 595 ! … … 597 597 ! 598 598 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 599 600 ! set transfer coefficients to default sea-ice values 601 Cd_atm(:,:) = Cd_ice 602 Ch_atm(:,:) = Cd_ice 603 Ce_atm(:,:) = Cd_ice 604 605 wndm_ice(:,:) = 0._wp !!gm brutal.... 606 607 ! ------------------------------------------------------------ ! 608 ! Wind module relative to the moving ice ( U10m - U_ice ) ! 609 ! ------------------------------------------------------------ ! 610 SELECT CASE( cp_ice_msh ) 611 CASE( 'I' ) ! B-grid ice dynamics : I-point (i.e. F-point with sea-ice indexation) 612 ! and scalar wind at T-point ( = | U10m - U_ice | ) (masked) 613 DO jj = 2, jpjm1 614 DO ji = 2, jpim1 ! B grid : NO vector opt 615 ! ... scalar wind at T-point (fld being at T-point) 616 zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.25 * ( u_ice(ji,jj+1) + u_ice(ji+1,jj+1) & 617 & + u_ice(ji,jj ) + u_ice(ji+1,jj ) ) 618 zwndj_t = sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.25 * ( v_ice(ji,jj+1) + v_ice(ji+1,jj+1) & 619 & + v_ice(ji,jj ) + v_ice(ji+1,jj ) ) 620 wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 621 END DO 622 END DO 623 CALL lbc_lnk( wndm_ice, 'T', 1. ) 624 ! 625 CASE( 'C' ) ! C-grid ice dynamics : U & V-points (same as ocean) 626 DO jj = 2, jpjm1 627 DO ji = fs_2, fs_jpim1 ! vect. opt. 628 zwndi_t = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) ) 629 zwndj_t = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji ,jj-1) + v_ice(ji,jj) ) ) 630 wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 631 END DO 632 END DO 633 CALL lbc_lnk( wndm_ice, 'T', 1. ) 634 ! 635 END SELECT 636 637 ! Make ice-atm. drag dependent on ice concentration 638 IF ( ln_Cd_L12 ) THEN ! calculate new drag from Lupkes(2012) equations 639 CALL Cdn10_Lupkes2012( Cd_atm ) 640 Ch_atm(:,:) = Cd_atm(:,:) ! momentum and heat transfer coef. are considered identical 641 ELSEIF( ln_Cd_L15 ) THEN ! calculate new drag from Lupkes(2015) equations 642 CALL Cdn10_Lupkes2015( Cd_atm, Ch_atm ) 607 643 ENDIF 608 #endif 644 645 !! CALL iom_put( "Cd_ice", Cd_atm) ! output value of pure ice-atm. transfer coef. 646 !! CALL iom_put( "Ch_ice", Ch_atm) ! output value of pure ice-atm. transfer coef. 609 647 610 648 ! local scalars ( place there for vector optimisation purposes) 611 649 ! Computing density of air! Way denser that 1.2 over sea-ice !!! 612 !!613 650 zrhoa (:,:) = rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) 614 651 … … 616 653 utau_ice (:,:) = 0._wp 617 654 vtau_ice (:,:) = 0._wp 618 wndm_ice (:,:) = 0._wp619 655 !!gm end 620 656 621 ! ------------------------------------------------------------ -----------------!622 ! Wind components and module relative to the moving ocean( U10m - U_ice ) !623 ! ------------------------------------------------------------ -----------------!657 ! ------------------------------------------------------------ ! 658 ! Wind stress relative to the moving ice ( U10m - U_ice ) ! 659 ! ------------------------------------------------------------ ! 624 660 SELECT CASE( cp_ice_msh ) 625 661 CASE( 'I' ) ! B-grid ice dynamics : I-point (i.e. F-point with sea-ice indexation) 626 ! and scalar wind at T-point ( = | U10m - U_ice | ) (masked)627 662 DO jj = 2, jpjm1 628 663 DO ji = 2, jpim1 ! B grid : NO vector opt … … 632 667 zwndj_f = 0.25 * ( sf(jp_wndj)%fnow(ji-1,jj ,1) + sf(jp_wndj)%fnow(ji ,jj ,1) & 633 668 & + 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 669 ! ... ice stress at I-point 670 zwnorm_f = zrhoa(ji,jj) * Cd_atm(ji,jj) * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f ) 636 671 utau_ice(ji,jj) = zwnorm_f * zwndi_f 637 672 vtau_ice(ji,jj) = zwnorm_f * zwndj_f 638 ! ... scalar wind at T-point (fld being at T-point)639 zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.25 * ( u_ice(ji,jj+1) + u_ice(ji+1,jj+1) &640 & + u_ice(ji,jj ) + u_ice(ji+1,jj ) )641 zwndj_t = sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.25 * ( v_ice(ji,jj+1) + v_ice(ji+1,jj+1) &642 & + v_ice(ji,jj ) + v_ice(ji+1,jj ) )643 wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)644 673 END DO 645 674 END DO 646 675 CALL lbc_lnk( utau_ice, 'I', -1. ) 647 676 CALL lbc_lnk( vtau_ice, 'I', -1. ) 648 CALL lbc_lnk( wndm_ice, 'T', 1. )649 677 ! 650 678 CASE( 'C' ) ! C-grid ice dynamics : U & V-points (same as ocean) 651 DO jj = 2, jpj652 DO ji = fs_2, jpi ! vect. opt.653 zwndi_t = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) )654 zwndj_t = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji ,jj-1) + v_ice(ji,jj) ) )655 wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)656 END DO657 END DO658 679 DO jj = 2, jpjm1 659 680 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) )&681 utau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji+1,jj ) + wndm_ice(ji,jj) ) & 661 682 & * ( 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) )&683 vtau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji,jj+1 ) + wndm_ice(ji,jj) ) & 663 684 & * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) ) 664 685 END DO … … 666 687 CALL lbc_lnk( utau_ice, 'U', -1. ) 667 688 CALL lbc_lnk( vtau_ice, 'V', -1. ) 668 CALL lbc_lnk( wndm_ice, 'T', 1. )669 689 ! 670 690 END SELECT … … 680 700 681 701 682 SUBROUTINE blk_ice_flx( ptsu, p alb )702 SUBROUTINE blk_ice_flx( ptsu, phs, phi, palb ) 683 703 !!--------------------------------------------------------------------- 684 704 !! *** ROUTINE blk_ice_flx *** … … 693 713 !!--------------------------------------------------------------------- 694 714 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: ptsu ! sea ice surface temperature 715 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phs ! snow thickness 716 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phi ! ice thickness 695 717 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: palb ! ice albedo (all skies) 696 718 !! … … 699 721 REAL(wp) :: zcoef_dqlw, zcoef_dqla ! - - 700 722 REAL(wp) :: zztmp, z1_lsub ! - - 723 REAL(wp) :: zfrqsr_tr_i ! sea ice shortwave fraction transmitted below through the ice 724 REAL(wp) :: zfr1, zfr2, zfac ! local variables 701 725 REAL(wp), DIMENSION(:,:,:), POINTER :: z_qlw ! long wave heat flux over ice 702 726 REAL(wp), DIMENSION(:,:,:), POINTER :: z_qsb ! sensible heat flux over ice … … 705 729 REAL(wp), DIMENSION(:,:) , POINTER :: zevap, zsnw ! evaporation and snw distribution after wind blowing (LIM3) 706 730 REAL(wp), DIMENSION(:,:) , POINTER :: zrhoa 707 REAL(wp), DIMENSION(:,:) , POINTER :: Cd ! transfer coefficient for momentum (tau) 731 708 732 !!--------------------------------------------------------------------- 709 733 ! … … 711 735 ! 712 736 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) 737 CALL wrk_alloc( jpi,jpj, zrhoa ) 738 ! 739 ! local scalars 727 740 zcoef_dqlw = 4.0 * 0.95 * Stef 728 741 zcoef_dqla = -Ls * 11637800. * (-5897.8) … … 752 765 ! ----------------------------! 753 766 754 ! ... turbulent heat fluxes 767 ! ... turbulent heat fluxes with Ch_atm recalculated in blk_ice_tau 755 768 ! 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))769 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 770 ! 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)) )771 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls * Ch_atm(ji,jj) * wndm_ice(ji,jj) * & 772 & ( 11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) ) 760 773 ! Latent heat sensitivity for ice (Dqla/Dt) 761 774 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))775 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 776 ELSE 764 777 dqla_ice(ji,jj,jl) = 0._wp … … 766 779 767 780 ! 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)781 z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj) 769 782 770 783 ! ----------------------------! … … 786 799 CALL iom_put( 'precip' , tprecip * 86400. ) ! Total precipitation 787 800 788 #if defined key_lim3789 801 CALL wrk_alloc( jpi,jpj, zevap, zsnw ) 790 802 … … 797 809 ! --- evaporation minus precipitation --- ! 798 810 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 )811 CALL ice_thd_snwblow( (1.-at_i_b(:,:)), zsnw ) ! snow distribution over ice after wind blowing 812 emp_oce(:,:) = ( 1._wp - at_i_b(:,:) ) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw ) 801 813 emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw 802 814 emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:) 803 815 804 816 ! --- heat flux associated with emp --- ! 805 qemp_oce(:,:) = - pfrld(:,:) * zevap(:,:) * sst_m(:,:) * rcp& ! evap at sst817 qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp & ! evap at sst 806 818 & + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp & ! liquid precip at Tair 807 819 & + sprecip(:,:) * ( 1._wp - zsnw ) * & ! solid precip at min(Tair,Tsnow) … … 811 823 812 824 ! --- 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 ) 825 qns_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) & 826 & + qemp_ice(:,:) + qemp_oce(:,:) 827 qsr_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 ) 815 828 816 829 ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- ! … … 824 837 825 838 CALL wrk_dealloc( jpi,jpj, zevap, zsnw ) 826 #endif 827 828 !-------------------------------------------------------------------- 829 ! FRACTIONs of net shortwave radiation which is not absorbed in the 830 ! thin surface layer and penetrates inside the ice cover 831 ! ( Maykut and Untersteiner, 1971 ; Ebert and Curry, 1993 ) 832 ! 833 fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 834 fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 835 ! 836 ! 839 840 ! --- absorbed and transmitted shortwave radiation (W/m2) --- ! 841 ! 842 ! former coding was 843 ! fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 844 ! fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 845 846 ! --- surface transmission parameter (i0, Grenfell Maykut 77) --- ! 847 zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) ! standard value 848 zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) ! zfr2 such that zfr1 + zfr2 to equal 1 849 850 qsr_ice_tr(:,:,:) = 0._wp 851 852 DO jl = 1, jpl 853 DO jj = 1, jpj 854 DO ji = 1, jpi 855 856 zfac = MAX( 0._wp , 1._wp - phi(ji,jj,jl) * 10._wp ) ! linear weighting factor: =0 for phi=0, =1 for phi = 0.1 857 zfrqsr_tr_i = zfr1 + zfac * zfr2 ! below 10 cm, linearly increase zfrqsr_tr_i until 1 at zero thickness 858 859 IF ( phs(ji,jj,jl) <= 0.0_wp ) THEN 860 zfrqsr_tr_i = zfr1 + zfac * zfr2 861 ELSE 862 zfrqsr_tr_i = 0._wp ! snow fully opaque 863 ENDIF 864 865 qsr_ice_tr(ji,jj,jl) = zfrqsr_tr_i * qsr_ice(ji,jj,jl) ! transmitted solar radiation 866 867 END DO 868 END DO 869 END DO 870 871 837 872 IF(ln_ctl) THEN 838 873 CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice : ', tab3d_2=z_qsb , clinfo2=' z_qsb : ', kdim=jpl) … … 846 881 CALL wrk_dealloc( jpi,jpj,jpl, z_qlw, z_qsb, z_dqlw, z_dqsb ) 847 882 CALL wrk_dealloc( jpi,jpj, zrhoa ) 848 CALL wrk_dealloc( jpi,jpj, Cd )849 883 ! 850 884 IF( nn_timing == 1 ) CALL timing_stop('blk_ice_flx') 851 885 852 886 END SUBROUTINE blk_ice_flx 887 888 889 890 SUBROUTINE blk_ice_qcn( k_monocat, ptsu, ptb, phs, phi ) 891 892 !!--------------------------------------------------------------------- 893 !! *** ROUTINE blk_ice_qcn *** 894 !! 895 !! ** Purpose : Compute surface temperature and snow/ice conduction flux 896 !! to force sea ice / snow thermodynamics 897 !! in the case JULES coupler is emulated 898 !! 899 !! ** Method : compute surface energy balance assuming neglecting heat storage 900 !! following the 0-layer Semtner (1976) approach 901 !! 902 !! ** Outputs : - ptsu : sea-ice / snow surface temperature (K) 903 !! - qcn_ice : surface inner conduction flux (W/m2) 904 !! 905 !!--------------------------------------------------------------------- 906 !! 907 INTEGER, INTENT(in) :: k_monocat ! single-category option 908 909 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ptsu ! sea ice / snow surface temperature 910 911 REAL(wp), DIMENSION(:,:) , INTENT(in) :: ptb ! sea ice base temperature 912 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phs ! snow thickness 913 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phi ! sea ice thickness 914 915 INTEGER :: ji, jj, jl ! dummy loop indices 916 INTEGER :: iter ! 917 REAL(wp) :: zfac, zfac2, zfac3 ! dummy factors 918 REAL(wp) :: zkeff_h, ztsu ! 919 REAL(wp) :: zqc, zqnet ! 920 REAL(wp) :: zhe, zqa0 ! 921 922 INTEGER , PARAMETER :: nit = 10 ! number of iterations 923 REAL(wp), PARAMETER :: zepsilon = 0.1_wp ! characteristic thickness for enhanced conduction 924 925 REAL(wp), DIMENSION(:,:,:), POINTER :: zgfac ! enhanced conduction factor 926 927 !!--------------------------------------------------------------------- 928 929 IF( nn_timing == 1 ) CALL timing_start('blk_ice_qcn') 930 ! 931 CALL wrk_alloc( jpi,jpj,jpl, zgfac ) 932 933 ! -------------------------------------! 934 ! I Enhanced conduction factor ! 935 ! -------------------------------------! 936 ! 937 ! Emulates the enhancement of conduction by unresolved thin ice (k_monocat = 1/3) 938 ! Fichefet and Morales Maqueda, JGR 1997 939 ! 940 zgfac(:,:,:) = 1._wp 941 942 SELECT CASE ( k_monocat ) 943 944 CASE ( 1 , 3 ) 945 946 zfac = 1._wp / ( rn_cnd_s + rcdic ) 947 zfac2 = EXP(1._wp) * 0.5_wp * zepsilon 948 zfac3 = 2._wp / zepsilon 949 950 DO jl = 1, jpl 951 DO jj = 1 , jpj 952 DO ji = 1, jpi 953 ! ! Effective thickness 954 zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcdic * phs(ji,jj,jl) ) * zfac 955 ! ! Enhanced conduction factor 956 IF( zhe >= zfac2 ) & 957 zgfac(ji,jj,jl) = MIN( 2._wp, ( 0.5_wp + 0.5 * LOG( zhe * zfac3 ) ) ) 958 END DO 959 END DO 960 END DO 961 962 END SELECT 963 964 ! -------------------------------------------------------------! 965 ! II Surface temperature and conduction flux ! 966 ! -------------------------------------------------------------! 967 968 zfac = rcdic * rn_cnd_s 969 ! ========================== ! 970 DO jl = 1, jpl ! Loop over ice categories ! 971 ! ! ========================== ! 972 DO jj = 1 , jpj 973 DO ji = 1, jpi 974 ! ! Effective conductivity of the snow-ice system divided by thickness 975 zkeff_h = zfac * zgfac(ji,jj,jl) / ( rcdic * phs(ji,jj,jl) + rn_cnd_s * phi(ji,jj,jl) ) 976 ! Store initial surface temperature 977 ztsu = ptsu(ji,jj,jl) 978 ! Net initial atmospheric heat flux 979 zqa0 = qsr_ice(ji,jj,jl) - qsr_ice_tr(ji,jj,jl) + qns_ice(ji,jj,jl) 980 981 DO iter = 1, nit ! --- Iteration loop 982 983 ! ! Conduction heat flux through snow-ice system (>0 downwards) 984 zqc = zkeff_h * ( ztsu - ptb(ji,jj) ) 985 ! ! Surface energy budget 986 zqnet = zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc 987 ! ! Temperature update 988 ztsu = ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h ) 989 990 END DO 991 992 ptsu(ji,jj,jl) = MIN( rt0, ztsu ) 993 994 qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) ) 995 996 END DO 997 END DO 998 999 END DO 1000 1001 CALL wrk_dealloc( jpi,jpj,jpl, zgfac ) 1002 1003 IF( nn_timing == 1 ) CALL timing_stop('blk_ice_qcn') 1004 1005 END SUBROUTINE blk_ice_qcn 853 1006 854 1007 #endif … … 973 1126 974 1127 #if defined key_lim3 1128 975 1129 SUBROUTINE Cdn10_Lupkes2012( Cd ) 976 1130 !!---------------------------------------------------------------------- … … 1022 1176 1023 1177 END SUBROUTINE Cdn10_Lupkes2012 1178 1179 1180 SUBROUTINE Cdn10_Lupkes2015( Cd, Ch ) 1181 !!---------------------------------------------------------------------- 1182 !! *** ROUTINE Cdn10_Lupkes2015 *** 1183 !! 1184 !! ** pUrpose : 1lternative turbulent transfert coefficients formulation 1185 !! between sea-ice and atmosphere with distinct momentum 1186 !! and heat coefficients depending on sea-ice concentration 1187 !! and atmospheric stability (no meltponds effect for now). 1188 !! 1189 !! ** Method : The parameterization is adapted from Lupkes et al. (2015) 1190 !! and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme, 1191 !! it considers specific skin and form drags (Andreas et al. 2010) 1192 !! to compute neutral transfert coefficients for both heat and 1193 !! momemtum fluxes. Atmospheric stability effect on transfert 1194 !! coefficient is also taken into account following Louis (1979). 1195 !! 1196 !! ** References : Lupkes et al. JGR 2015 (theory) 1197 !! Lupkes et al. ECHAM6 documentation 2015 (implementation) 1198 !! 1199 !!---------------------------------------------------------------------- 1200 ! 1201 REAL(wp), DIMENSION(:,:), INTENT(inout) :: Cd 1202 REAL(wp), DIMENSION(:,:), INTENT(inout) :: Ch 1203 REAL(wp), DIMENSION(jpi,jpj) :: zst, zqo_sat, zqi_sat 1204 ! 1205 ! ECHAM6 constants 1206 REAL(wp), PARAMETER :: z0_skin_ice = 0.69e-3_wp ! Eq. 43 [m] 1207 REAL(wp), PARAMETER :: z0_form_ice = 0.57e-3_wp ! Eq. 42 [m] 1208 REAL(wp), PARAMETER :: z0_ice = 1.00e-3_wp ! Eq. 15 [m] 1209 REAL(wp), PARAMETER :: zce10 = 2.80e-3_wp ! Eq. 41 1210 REAL(wp), PARAMETER :: zbeta = 1.1_wp ! Eq. 41 1211 REAL(wp), PARAMETER :: zc = 5._wp ! Eq. 13 1212 REAL(wp), PARAMETER :: zc2 = zc * zc 1213 REAL(wp), PARAMETER :: zam = 2. * zc ! Eq. 14 1214 REAL(wp), PARAMETER :: zah = 3. * zc ! Eq. 30 1215 REAL(wp), PARAMETER :: z1_alpha = 1._wp / 0.2_wp ! Eq. 51 1216 REAL(wp), PARAMETER :: z1_alphaf = z1_alpha ! Eq. 56 1217 REAL(wp), PARAMETER :: zbetah = 1.e-3_wp ! Eq. 26 1218 REAL(wp), PARAMETER :: zgamma = 1.25_wp ! Eq. 26 1219 REAL(wp), PARAMETER :: z1_gamma = 1._wp / zgamma 1220 REAL(wp), PARAMETER :: r1_3 = 1._wp / 3._wp 1221 ! 1222 INTEGER :: ji, jj ! dummy loop indices 1223 REAL(wp) :: zthetav_os, zthetav_is, zthetav_zu 1224 REAL(wp) :: zrib_o, zrib_i 1225 REAL(wp) :: zCdn_skin_ice, zCdn_form_ice, zCdn_ice 1226 REAL(wp) :: zChn_skin_ice, zChn_form_ice 1227 REAL(wp) :: z0w, z0i, zfmi, zfmw, zfhi, zfhw 1228 REAL(wp) :: zCdn_form_tmp 1229 !!---------------------------------------------------------------------- 1230 1231 ! Momentum Neutral Transfert Coefficients (should be a constant) 1232 zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2 ! Eq. 40 1233 zCdn_skin_ice = ( vkarmn / LOG( rn_zu / z0_skin_ice + 1._wp ) )**2 ! Eq. 7 1234 zCdn_ice = zCdn_skin_ice ! Eq. 7 (cf Lupkes email for details) 1235 !zCdn_ice = 1.89e-3 ! old ECHAM5 value (cf Eq. 32) 1236 1237 ! Heat Neutral Transfert Coefficients 1238 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) 1239 1240 ! Atmospheric and Surface Variables 1241 zst(:,:) = sst_m(:,:) + rt0 ! convert SST from Celcius to Kelvin 1242 zqo_sat(:,:) = 0.98_wp * q_sat( zst(:,:) , sf(jp_slp)%fnow(:,:,1) ) ! saturation humidity over ocean [kg/kg] 1243 zqi_sat(:,:) = 0.98_wp * q_sat( tm_su(:,:), sf(jp_slp)%fnow(:,:,1) ) ! saturation humidity over ice [kg/kg] 1244 ! 1245 DO jj = 2, jpjm1 ! reduced loop is necessary for reproducibility 1246 DO ji = fs_2, fs_jpim1 1247 ! Virtual potential temperature [K] 1248 zthetav_os = zst(ji,jj) * ( 1._wp + rctv0 * zqo_sat(ji,jj) ) ! over ocean 1249 zthetav_is = tm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) ) ! ocean ice 1250 zthetav_zu = t_zu (ji,jj) * ( 1._wp + rctv0 * q_zu(ji,jj) ) ! at zu 1251 1252 ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead) 1253 zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj) )**2 ! over ocean 1254 zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2 ! over ice 1255 1256 ! Momentum and Heat Neutral Transfert Coefficients 1257 zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta ! Eq. 40 1258 zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) ) ! Eq. 53 1259 1260 ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead) 1261 z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water 1262 z0i = z0_skin_ice ! over ice (cf Lupkes email for details) 1263 IF( zrib_o <= 0._wp ) THEN 1264 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 1265 zfhw = ( 1._wp + ( zbetah * ( zthetav_os - zthetav_zu )**r1_3 / ( Chn_oce(ji,jj) * MAX(0.01, wndm(ji,jj)) ) & ! Eq. 26 1266 & )**zgamma )**z1_gamma 1267 ELSE 1268 zfmw = 1._wp / ( 1._wp + zam * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 12 1269 zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 28 1270 ENDIF 1271 1272 IF( zrib_i <= 0._wp ) THEN 1273 zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 9 1274 zfhi = 1._wp - zah * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 25 1275 ELSE 1276 zfmi = 1._wp / ( 1._wp + zam * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 11 1277 zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 27 1278 ENDIF 1279 1280 ! Momentum Transfert Coefficients (Eq. 38) 1281 Cd(ji,jj) = zCdn_skin_ice * zfmi + & 1282 & 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) ) 1283 1284 ! Heat Transfert Coefficients (Eq. 49) 1285 Ch(ji,jj) = zChn_skin_ice * zfhi + & 1286 & 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) ) 1287 ! 1288 END DO 1289 END DO 1290 CALL lbc_lnk_multi( Cd, 'T', 1., Ch, 'T', 1. ) 1291 ! 1292 END SUBROUTINE Cdn10_Lupkes2015 1293 1024 1294 #endif 1025 1026 1295 1027 1296 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
