MODULE sbcblk !!====================================================================== !! *** MODULE sbcblk *** !! Ocean forcing: momentum, heat and freshwater flux formulation !! Aerodynamic Bulk Formulas !! SUCCESSOR OF "sbcblk_core" !!===================================================================== !! History : 1.0 ! 2004-08 (U. Schweckendiek) Original CORE code !! 2.0 ! 2005-04 (L. Brodeau, A.M. Treguier) improved CORE bulk and its user interface !! 3.0 ! 2006-06 (G. Madec) sbc rewritting !! - ! 2006-12 (L. Brodeau) Original code for turb_core !! 3.2 ! 2009-04 (B. Lemaire) Introduce iom_put !! 3.3 ! 2010-10 (S. Masson) add diurnal cycle !! 3.4 ! 2011-11 (C. Harris) Fill arrays required by CICE !! 3.7 ! 2014-06 (L. Brodeau) simplification and optimization of CORE bulk !! 4.0 ! 2016-06 (L. Brodeau) sbcblk_core becomes sbcblk and is not restricted to the CORE algorithm anymore !! ! ==> based on AeroBulk (http://aerobulk.sourceforge.net/) !! 4.0 ! 2016-10 (G. Madec) introduce a sbc_blk_init routine !! 4.0 ! 2016-10 (M. Vancoppenolle) Introduce conduction flux emulator (M. Vancoppenolle) !! 4.0 ! 2019-03 (F. Lemarié & G. Samson) add ABL compatibility (ln_abl=TRUE) !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! sbc_blk_init : initialisation of the chosen bulk formulation as ocean surface boundary condition !! sbc_blk : bulk formulation as ocean surface boundary condition !! blk_oce_1 : computes pieces of momentum, heat and freshwater fluxes over ocean for ABL model (ln_abl=TRUE) !! blk_oce_2 : finalizes momentum, heat and freshwater fluxes computation over ocean after the ABL step (ln_abl=TRUE) !! rho_air : density of (moist) air (depends on T_air, q_air and SLP !! cp_air : specific heat of (moist) air (depends spec. hum. q_air) !! q_sat : saturation humidity as a function of SLP and temperature !! L_vap : latent heat of vaporization of water as a function of temperature !! sea-ice case only : !! blk_ice_1 : provide the air-ice stress !! blk_ice_2 : provide the heat and mass fluxes at air-ice interface !! blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux) !! Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag !! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag !!---------------------------------------------------------------------- USE oce ! ocean dynamics and tracers USE dom_oce ! ocean space and time domain USE phycst ! physical constants USE fldread ! read input fields USE sbc_oce ! Surface boundary condition: ocean fields USE cyclone ! Cyclone 10m wind form trac of cyclone centres USE sbcdcy ! surface boundary condition: diurnal cycle USE sbcwave , ONLY : cdn_wave ! wave module USE sbc_ice ! Surface boundary condition: ice fields USE lib_fortran ! to use key_nosignedzero #if defined key_si3 USE ice , ONLY : jpl, a_i_b, at_i_b, rn_cnd_s, hfx_err_dif USE icethd_dh ! for CALL ice_thd_snwblow #endif USE sbcblk_algo_ncar ! => turb_ncar : NCAR - CORE (Large & Yeager, 2009) USE sbcblk_algo_coare ! => turb_coare : COAREv3.0 (Fairall et al. 2003) USE sbcblk_algo_coare3p5 ! => turb_coare3p5 : COAREv3.5 (Edson et al. 2013) USE sbcblk_algo_ecmwf ! => turb_ecmwf : ECMWF (IFS cycle 31) ! USE iom ! I/O manager library USE in_out_manager ! I/O manager USE lib_mpp ! distribued memory computing library USE lbclnk ! ocean lateral boundary conditions (or mpp link) USE prtctl ! Print control IMPLICIT NONE PRIVATE PUBLIC sbc_blk_init ! called in sbcmod PUBLIC sbc_blk ! called in sbcmod PUBLIC blk_oce_1 ! called in sbcabl PUBLIC blk_oce_2 ! called in sbcabl PUBLIC rho_air ! called in ablmod PUBLIC cp_air ! called in ablmod #if defined key_si3 PUBLIC blk_ice_1 ! routine called in icesbc PUBLIC blk_ice_2 ! routine called in icesbc PUBLIC blk_ice_qcn ! routine called in icesbc #endif INTERFACE cp_air MODULE PROCEDURE cp_air_0d, cp_air_2d END INTERFACE ! !!* Namelist namsbc_blk : bulk parameters LOGICAL :: ln_NCAR ! "NCAR" algorithm (Large and Yeager 2008) LOGICAL :: ln_COARE_3p0 ! "COARE 3.0" algorithm (Fairall et al. 2003) LOGICAL :: ln_COARE_3p5 ! "COARE 3.5" algorithm (Edson et al. 2013) LOGICAL :: ln_ECMWF ! "ECMWF" algorithm (IFS cycle 31) ! ! REAL(wp) :: rn_pfac ! multiplication factor for precipitation REAL(wp), PUBLIC :: rn_efac !: multiplication factor for evaporation REAL(wp), PUBLIC :: rn_vfac !: multiplication factor for ice/ocean velocity in the calculation of wind stress REAL(wp) :: rn_zqt ! z(q,t) : height of humidity and temperature measurements REAL(wp) :: rn_zu ! z(u) : height of wind measurements ! ! LOGICAL :: ln_Cd_L12 ! ice-atm drag = F( ice concentration ) (Lupkes et al. JGR2012) LOGICAL :: ln_Cd_L15 ! ice-atm drag = F( ice concentration, atmospheric stability ) (Lupkes et al. JGR2015) INTEGER :: nblk ! choice of the bulk algorithm ! ! associated indices: INTEGER, PARAMETER :: np_NCAR = 1 ! "NCAR" algorithm (Large and Yeager 2008) INTEGER, PARAMETER :: np_COARE_3p0 = 2 ! "COARE 3.0" algorithm (Fairall et al. 2003) INTEGER, PARAMETER :: np_COARE_3p5 = 3 ! "COARE 3.5" algorithm (Edson et al. 2013) INTEGER, PARAMETER :: np_ECMWF = 4 ! "ECMWF" algorithm (IFS cycle 31) ! !!! air parameters REAL(wp) , PARAMETER :: Cp_dry = 1005.0 ! Specic heat of dry air, constant pressure [J/K/kg] REAL(wp) , PARAMETER :: Cp_vap = 1860.0 ! Specic heat of water vapor, constant pressure [J/K/kg] REAL(wp), PUBLIC, PARAMETER :: R_dry = 287.05_wp !: Specific gas constant for dry air [J/K/kg] REAL(wp) , PARAMETER :: R_vap = 461.495_wp ! Specific gas constant for water vapor [J/K/kg] REAL(wp) , PARAMETER :: reps0 = R_dry/R_vap ! ratio of gas constant for dry air and water vapor => ~ 0.622 REAL(wp), PUBLIC, PARAMETER :: rctv0 = R_vap/R_dry !: for virtual temperature (== (1-eps)/eps) => ~ 0.608 ! !!! Bulk parameters REAL(wp) , PARAMETER :: cpa = 1000.5 ! specific heat of air (only used for ice fluxes now...) REAL(wp) , PARAMETER :: Ls = 2.839e6 ! latent heat of sublimation REAL(wp) , PARAMETER :: Stef = 5.67e-8 ! Stefan Boltzmann constant REAL(wp) , PARAMETER :: rcdice = 1.4e-3 ! transfer coefficient over ice REAL(wp) , PARAMETER :: albo = 0.066 ! ocean albedo assumed to be constant ! REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: Cd_ice , Ch_ice , Ce_ice ! transfert coefficients over ice REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: Cdn_oce, Chn_oce, Cen_oce ! neutral coeffs over ocean (L15 bulk scheme) REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: t_zu, q_zu ! air temp. and spec. hum. at wind speed height (L15 bulk scheme) !INTEGER , PUBLIC, PARAMETER :: jpfld =11 !: maximum number of files to read INTEGER , PUBLIC :: jpfld !: maximum number of files to read INTEGER , PUBLIC, PARAMETER :: jp_wndi = 1 !: index of 10m wind velocity (i-component) (m/s) at T-point INTEGER , PUBLIC, PARAMETER :: jp_wndj = 2 !: index of 10m wind velocity (j-component) (m/s) at T-point INTEGER , PUBLIC, PARAMETER :: jp_tair = 3 !: index of 10m air temperature (Kelvin) INTEGER , PUBLIC, PARAMETER :: jp_humi = 4 !: index of specific humidity ( % ) INTEGER , PUBLIC, PARAMETER :: jp_qsr = 5 !: index of solar heat (W/m2) INTEGER , PUBLIC, PARAMETER :: jp_qlw = 6 !: index of Long wave (W/m2) INTEGER , PUBLIC, PARAMETER :: jp_prec = 7 !: index of total precipitation (rain+snow) (Kg/m2/s) INTEGER , PUBLIC, PARAMETER :: jp_snow = 8 !: index of snow (solid prcipitation) (kg/m2/s) INTEGER , PUBLIC, PARAMETER :: jp_slp = 9 !: index of sea level pressure (Pa) INTEGER , PUBLIC, PARAMETER :: jp_hpgi =10 !: index of ABL geostrophic wind or hpg (i-component) (m/s) at T-point INTEGER , PUBLIC, PARAMETER :: jp_hpgj =11 !: index of ABL geostrophic wind or hpg (j-component) (m/s) at T-point TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) :: sf ! structure of input atmospheric fields (file informations, fields read) !! * Substitutions # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/OCE 4.0 , NEMO Consortium (2018) !! $Id$ !! Software governed by the CeCILL license (see ./LICENSE) !!---------------------------------------------------------------------- CONTAINS INTEGER FUNCTION sbc_blk_alloc() !!------------------------------------------------------------------- !! *** ROUTINE sbc_blk_alloc *** !!------------------------------------------------------------------- ALLOCATE( t_zu(jpi,jpj), q_zu(jpi,jpj), & & Cdn_oce(jpi,jpj), Chn_oce(jpi,jpj), Cen_oce(jpi,jpj), & & Cd_ice (jpi,jpj), Ch_ice (jpi,jpj), Ce_ice (jpi,jpj), STAT=sbc_blk_alloc ) ! CALL mpp_sum ( 'sbcblk', sbc_blk_alloc ) IF( sbc_blk_alloc /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk_alloc: failed to allocate arrays' ) END FUNCTION sbc_blk_alloc SUBROUTINE sbc_blk_init !!--------------------------------------------------------------------- !! *** ROUTINE sbc_blk_init *** !! !! ** Purpose : choose and initialize a bulk formulae formulation !! !! ** Method : !! !!---------------------------------------------------------------------- INTEGER :: jfpr ! dummy loop indice and argument INTEGER :: ios, ierror, ioptio ! Local integer !! CHARACTER(len=100) :: cn_dir ! Root directory for location of atmospheric forcing files !TYPE(FLD_N), DIMENSION(jpfld) :: slf_i ! array of namelist informations on the fields to read TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) :: slf_i ! array of namelist informations on the fields to read TYPE(FLD_N) :: sn_wndi, sn_wndj, sn_humi, sn_qsr ! informations about the fields to be read TYPE(FLD_N) :: sn_qlw , sn_tair, sn_prec, sn_snow ! " " TYPE(FLD_N) :: sn_slp , sn_hpgi, sn_hpgj ! " " NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw , & ! input fields & sn_tair, sn_prec, sn_snow, sn_slp, sn_hpgi, sn_hpgj, & & ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF, & ! bulk algorithm & cn_dir , rn_zqt, rn_zu, & & rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15 !!--------------------------------------------------------------------- ! ! ! allocate sbc_blk_core array IF( sbc_blk_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' ) ! ! !** read bulk namelist REWIND( numnam_ref ) !* Namelist namsbc_blk in reference namelist : bulk parameters READ ( numnam_ref, namsbc_blk, IOSTAT = ios, ERR = 901) 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_blk in reference namelist' ) ! REWIND( numnam_cfg ) !* Namelist namsbc_blk in configuration namelist : bulk parameters READ ( numnam_cfg, namsbc_blk, IOSTAT = ios, ERR = 902 ) 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namsbc_blk in configuration namelist' ) ! IF(lwm) WRITE( numond, namsbc_blk ) ! ! !** initialization of the chosen bulk formulae (+ check) ! !* select the bulk chosen in the namelist and check the choice ioptio = 0 IF( ln_NCAR ) THEN ; nblk = np_NCAR ; ioptio = ioptio + 1 ; ENDIF IF( ln_COARE_3p0 ) THEN ; nblk = np_COARE_3p0 ; ioptio = ioptio + 1 ; ENDIF IF( ln_COARE_3p5 ) THEN ; nblk = np_COARE_3p5 ; ioptio = ioptio + 1 ; ENDIF IF( ln_ECMWF ) THEN ; nblk = np_ECMWF ; ioptio = ioptio + 1 ; ENDIF ! IF( ioptio /= 1 ) CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' ) ! IF( ln_dm2dc ) THEN !* check: diurnal cycle on Qsr IF( sn_qsr%freqh /= 24. ) CALL ctl_stop( 'sbc_blk_init: ln_dm2dc=T only with daily short-wave input' ) IF( sn_qsr%ln_tint ) THEN CALL ctl_warn( 'sbc_blk_init: ln_dm2dc=T daily qsr time interpolation done by sbcdcy module', & & ' ==> We force time interpolation = .false. for qsr' ) sn_qsr%ln_tint = .false. ENDIF ENDIF ! !* set the bulk structure ! !- store namelist information in an array IF( ln_blk ) jpfld = 9 IF( ln_abl ) jpfld = 11 ALLOCATE( slf_i(jpfld) ) ! slf_i(jp_wndi) = sn_wndi ; slf_i(jp_wndj) = sn_wndj slf_i(jp_qsr ) = sn_qsr ; slf_i(jp_qlw ) = sn_qlw slf_i(jp_tair) = sn_tair ; slf_i(jp_humi) = sn_humi slf_i(jp_prec) = sn_prec ; slf_i(jp_snow) = sn_snow slf_i(jp_slp ) = sn_slp IF( ln_abl ) THEN slf_i(jp_hpgi) = sn_hpgi ; slf_i(jp_hpgj) = sn_hpgj END IF ! ! !- allocate the bulk structure ALLOCATE( sf(jpfld), STAT=ierror ) IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_blk_init: unable to allocate sf structure' ) ! !- fill the bulk structure with namelist informations CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' ) ! DO jfpr= 1, jpfld ! IF( TRIM(sf(jfpr)%clrootname) == 'NOT USED' ) THEN !-- not used field --! (only now allocated and set to zero) ALLOCATE( sf(jfpr)%fnow(jpi,jpj,1) ) sf(jfpr)%fnow(:,:,1) = 0._wp ELSE !-- used field --! IF( ln_abl .AND. & & ( jfpr == jp_wndi .OR. jfpr == jp_wndj .OR. jfpr == jp_humi .OR. & & jfpr == jp_hpgi .OR. jfpr == jp_hpgj .OR. jfpr == jp_tair ) ) THEN ! ABL: some fields are 3D input ALLOCATE( sf(jfpr)%fnow(jpi,jpj,jpka) ) IF( slf_i(jfpr)%ln_tint ) ALLOCATE( sf(jfpr)%fdta(jpi,jpj,jpka,2) ) ELSE ! others or Bulk fields are 2D fiels ALLOCATE( sf(jfpr)%fnow(jpi,jpj,1) ) IF( slf_i(jfpr)%ln_tint ) ALLOCATE( sf(jfpr)%fdta(jpi,jpj,1,2) ) ENDIF ! IF( slf_i(jfpr)%freqh > 0. .AND. MOD( NINT(3600. * slf_i(jfpr)%freqh), nn_fsbc * NINT(rdt) ) /= 0 ) & & CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rdt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.', & & ' This is not ideal. You should consider changing either rdt or nn_fsbc value...' ) ENDIF END DO ! IF( ln_wave ) THEN ! surface waves IF( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) ) & ! Activated wave module but neither drag nor stokes drift activated & CALL ctl_stop( 'sbc_blk_init: Ask for wave coupling but ln_cdgw=ln_sdw=ln_tauwoc=ln_stcor=F' ) IF( ln_cdgw .AND. .NOT.ln_NCAR ) & ! drag coefficient read from wave model only with NCAR bulk formulae & CALL ctl_stop( 'sbc_blk_init: drag coefficient read from wave model need NCAR bulk formulae') IF( ln_stcor .AND. .NOT.ln_sdw ) & CALL ctl_stop( 'sbc_blk_init: Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T') ELSE IF( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) THEN CALL ctl_warn( 'sbc_blk_init: ln_wave=F, set all wave-related namelist parameter to FALSE') ln_cdgw =.FALSE. ; ln_sdw =.FALSE. ; ln_tauwoc =.FALSE. ; ln_stcor =.FALSE. ENDIF ENDIF ! IF( ln_abl ) THEN ! ABL: read 3D fields for wind, temperature, humidity and pressure gradient rn_zqt = ght_abl(2) ! set the bulk altitude to ABL first level rn_zu = ght_abl(2) IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) ' ABL formulation: overwrite rn_zqt & rn_zu with ABL first level altitude' ENDIF ! ! set transfer coefficients to default sea-ice values Cd_ice(:,:) = rcdice Ch_ice(:,:) = rcdice Ce_ice(:,:) = rcdice ! IF(lwp) THEN !** Control print ! WRITE(numout,*) !* namelist WRITE(numout,*) ' Namelist namsbc_blk (other than data information):' WRITE(numout,*) ' "NCAR" algorithm (Large and Yeager 2008) ln_NCAR = ', ln_NCAR WRITE(numout,*) ' "COARE 3.0" algorithm (Fairall et al. 2003) ln_COARE_3p0 = ', ln_COARE_3p0 WRITE(numout,*) ' "COARE 3.5" algorithm (Edson et al. 2013) ln_COARE_3p5 = ', ln_COARE_3p0 WRITE(numout,*) ' "ECMWF" algorithm (IFS cycle 31) ln_ECMWF = ', ln_ECMWF WRITE(numout,*) ' Air temperature and humidity reference height (m) rn_zqt = ', rn_zqt WRITE(numout,*) ' Wind vector reference height (m) rn_zu = ', rn_zu WRITE(numout,*) ' factor applied on precipitation (total & snow) rn_pfac = ', rn_pfac WRITE(numout,*) ' factor applied on evaporation rn_efac = ', rn_efac WRITE(numout,*) ' factor applied on ocean/ice velocity rn_vfac = ', rn_vfac WRITE(numout,*) ' (form absolute (=0) to relative winds(=1))' WRITE(numout,*) ' use ice-atm drag from Lupkes2012 ln_Cd_L12 = ', ln_Cd_L12 WRITE(numout,*) ' use ice-atm drag from Lupkes2015 ln_Cd_L15 = ', ln_Cd_L15 ! WRITE(numout,*) SELECT CASE( nblk ) !* Print the choice of bulk algorithm CASE( np_NCAR ) ; WRITE(numout,*) ' ==>>> "NCAR" algorithm (Large and Yeager 2008)' CASE( np_COARE_3p0 ) ; WRITE(numout,*) ' ==>>> "COARE 3.0" algorithm (Fairall et al. 2003)' CASE( np_COARE_3p5 ) ; WRITE(numout,*) ' ==>>> "COARE 3.5" algorithm (Edson et al. 2013)' CASE( np_ECMWF ) ; WRITE(numout,*) ' ==>>> "ECMWF" algorithm (IFS cycle 31)' END SELECT ! ENDIF ! END SUBROUTINE sbc_blk_init SUBROUTINE sbc_blk( kt ) !!--------------------------------------------------------------------- !! *** ROUTINE sbc_blk *** !! !! ** Purpose : provide at each time step the surface ocean fluxes !! (momentum, heat, freshwater and runoff) !! !! ** Method : (1) READ each fluxes in NetCDF files: !! the 10m wind velocity (i-component) (m/s) at T-point !! the 10m wind velocity (j-component) (m/s) at T-point !! the 10m or 2m specific humidity (kg/kg) !! the solar heat (W/m2) !! the Long wave (W/m2) !! the 10m or 2m air temperature (Kelvin) !! the total precipitation (rain+snow) (Kg/m2/s) !! the snow (solid precipitation) (kg/m2/s) !! ABL dynamical forcing (i/j-components of either hpg or geostrophic winds) !! (2) CALL blk_oce_1 and blk_oce_2 !! !! C A U T I O N : never mask the surface stress fields !! the stress is assumed to be in the (i,j) mesh referential !! !! ** Action : defined at each time-step at the air-sea interface !! - utau, vtau i- and j-component of the wind stress !! - taum wind stress module at T-point !! - wndm wind speed module at T-point over free ocean or leads in presence of sea-ice !! - qns, qsr non-solar and solar heat fluxes !! - emp upward mass flux (evapo. - precip.) !! - sfx salt flux due to freezing/melting (non-zero only if ice is present) !! !! ** References : Large & Yeager, 2004 / Large & Yeager, 2008 !! Brodeau et al. Ocean Modelling 2010 !!---------------------------------------------------------------------- INTEGER, INTENT(in) :: kt ! ocean time step !! REAL(wp), DIMENSION(jpi,jpj) :: zssq, zcd_du, zsen, zevp !!--------------------------------------------------------------------- ! CALL fld_read( kt, nn_fsbc, sf ) ! input fields provided at the current time-step ! ! ! compute the surface ocean fluxes using bulk formulea IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN ! CALL blk_oce_1( kt, sf(jp_wndi)%fnow(:,:,1), sf(jp_wndj)%fnow(:,:,1), & ! <<= in & sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), & ! <<= in & sf(jp_slp )%fnow(:,:,1), sst_m, ssu_m, ssv_m, & ! <<= in & zssq, zcd_du, zsen, zevp ) ! =>> out CALL blk_oce_2( sf(jp_tair)%fnow(:,:,1), sf(jp_qsr )%fnow(:,:,1), & ! <<= in & sf(jp_qlw )%fnow(:,:,1), sf(jp_prec)%fnow(:,:,1), & ! <<= in & sf(jp_snow)%fnow(:,:,1), sst_m, & ! <<= in & zsen, zevp ) ! <=> in out ENDIF #if defined key_cice IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN qlw_ice(:,:,1) = sf(jp_qlw )%fnow(:,:,1) IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) ELSE ; qsr_ice(:,:,1) = sf(jp_qsr)%fnow(:,:,1) ENDIF tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) qatm_ice(:,:) = sf(jp_humi)%fnow(:,:,1) tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac wndi_ice(:,:) = sf(jp_wndi)%fnow(:,:,1) wndj_ice(:,:) = sf(jp_wndj)%fnow(:,:,1) ENDIF #endif ! END SUBROUTINE sbc_blk SUBROUTINE blk_oce_1( kt, pwndi, pwndj, ptair, phumi, & ! inp & pslp, pst , pu , pv, & ! inp & pssq, pcd_du, psen, pevp ) ! out !!--------------------------------------------------------------------- !! *** ROUTINE blk_oce_1 *** !! !! ** Purpose : Computes Cd x |U|, Ch x |U|, Ce x |U| for ABL integration !! !! ** Method : bulk formulae using atmospheric !! fields from the ABL model at previous time-step !! !! ** Outputs : - pssq : surface humidity used to compute latent heat flux (kg/kg) !! - pcd_du : Cd x |dU| at T-points (m/s) !! - psen : Ch x |dU| at T-points (m/s) !! - pevp : Ce x |dU| at T-points (m/s) !!--------------------------------------------------------------------- INTEGER , INTENT(in ) :: kt ! time step index REAL(wp), INTENT(in ), DIMENSION(:,:) :: pwndi ! atmospheric wind at U-point [m/s] REAL(wp), INTENT(in ), DIMENSION(:,:) :: pwndj ! atmospheric wind at V-point [m/s] REAL(wp), INTENT(in ), DIMENSION(:,:) :: phumi ! specific humidity at T-points [kg/kg] REAL(wp), INTENT(in ), DIMENSION(:,:) :: ptair ! potential temperature at T-points [Kelvin] REAL(wp), INTENT(in ), DIMENSION(:,:) :: pslp ! sea-level pressure [Pa] REAL(wp), INTENT(in ), DIMENSION(:,:) :: pst ! surface temperature [Celcius] REAL(wp), INTENT(in ), DIMENSION(:,:) :: pu ! surface current at U-point (i-component) [m/s] REAL(wp), INTENT(in ), DIMENSION(:,:) :: pv ! surface current at V-point (j-component) [m/s] REAL(wp), INTENT( out), DIMENSION(:,:) :: pssq ! specific humidity at pst [kg/kg] REAL(wp), INTENT( out), DIMENSION(:,:) :: pcd_du ! Cd x |dU| at T-points [m/s] REAL(wp), INTENT( out), DIMENSION(:,:) :: psen ! Ch x |dU| at T-points [m/s] REAL(wp), INTENT( out), DIMENSION(:,:) :: pevp ! Ce x |dU| at T-points [m/s] ! INTEGER :: ji, jj ! dummy loop indices REAL(wp) :: zztmp ! local variable REAL(wp), DIMENSION(jpi,jpj) :: zwnd_i, zwnd_j ! wind speed components at T-point REAL(wp), DIMENSION(jpi,jpj) :: zst ! surface temperature in Kelvin REAL(wp), DIMENSION(jpi,jpj) :: zU_zu ! bulk wind speed at height zu [m/s] REAL(wp), DIMENSION(jpi,jpj) :: ztpot ! potential temperature of air at z=rn_zqt [K] REAL(wp), DIMENSION(jpi,jpj) :: zrhoa ! density of air [kg/m^3] REAL(wp), DIMENSION(jpi,jpj) :: zcd_oce ! momentum transfert coefficient over ocean REAL(wp), DIMENSION(jpi,jpj) :: zch_oce ! sensible heat transfert coefficient over ocean REAL(wp), DIMENSION(jpi,jpj) :: zce_oce ! latent heat transfert coefficient over ocean !!--------------------------------------------------------------------- ! ! local scalars ( place there for vector optimisation purposes) zst(:,:) = pst(:,:) + rt0 ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K) ! ----------------------------------------------------------------------------- ! ! 0 Wind components and module at T-point relative to the moving ocean ! ! ----------------------------------------------------------------------------- ! ! ... components ( U10m - U_oce ) at T-point (unmasked) #if defined key_cyclone zwnd_i(:,:) = 0._wp zwnd_j(:,:) = 0._wp CALL wnd_cyc( kt, zwnd_i, zwnd_j ) ! add analytical tropical cyclone (Vincent et al. JGR 2012) DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vect. opt. pwndi(ji,jj) = pwndi(ji,jj) + zwnd_i(ji,jj) pwndj(ji,jj) = pwndj(ji,jj) + zwnd_j(ji,jj) END DO END DO #endif DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vect. opt. zwnd_i(ji,jj) = ( pwndi(ji,jj) - rn_vfac * 0.5 * ( pu(ji-1,jj ) + pu(ji,jj) ) ) zwnd_j(ji,jj) = ( pwndj(ji,jj) - rn_vfac * 0.5 * ( pv(ji ,jj-1) + pv(ji,jj) ) ) END DO END DO CALL lbc_lnk_multi( 'sbcblk', zwnd_i, 'T', -1., zwnd_j, 'T', -1. ) ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked) wndm(:,:) = SQRT( zwnd_i(:,:) * zwnd_i(:,:) & & + zwnd_j(:,:) * zwnd_j(:,:) ) * tmask(:,:,1) ! ----------------------------------------------------------------------------- ! ! I Turbulent FLUXES ! ! ----------------------------------------------------------------------------- ! ! ... specific humidity at SST and IST tmask( pssq(:,:) = 0.98 * q_sat( zst(:,:), pslp(:,:) ) ! IF( ln_abl ) THEN ztpot = ptair(:,:) ELSE ! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate ! (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2 ! (since reanalysis products provide T at z, not theta !) ztpot = ptair(:,:) + gamma_moist( ptair(:,:), phumi(:,:) ) * rn_zqt ENDIF SELECT CASE( nblk ) !== transfer coefficients ==! Cd, Ch, Ce at T-point ! CASE( np_NCAR ) ; CALL turb_ncar ( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm, & ! NCAR-COREv2 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, Cdn_oce, Chn_oce, Cen_oce ) CASE( np_COARE_3p0 ) ; CALL turb_coare ( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm, & ! COARE v3.0 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, Cdn_oce, Chn_oce, Cen_oce ) CASE( np_COARE_3p5 ) ; CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm, & ! COARE v3.5 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, Cdn_oce, Chn_oce, Cen_oce ) CASE( np_ECMWF ) ; CALL turb_ecmwf ( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm, & ! ECMWF & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, Cdn_oce, Chn_oce, Cen_oce ) CASE DEFAULT CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) END SELECT IF ( ln_abl ) THEN !== ABL formulation ==! multiplication by rho_air and turbulent fluxes computation done in ablstp !! FL do we need this multiplication by tmask ... ??? DO jj = 1, jpj DO ji = 1, jpi zztmp = zU_zu(ji,jj) !* tmask(ji,jj,1) wndm(ji,jj) = zztmp ! Store zU_zu in wndm to compute ustar2 in ablmod pcd_du(ji,jj) = zztmp * zcd_oce(ji,jj) psen(ji,jj) = zztmp * zch_oce(ji,jj) pevp(ji,jj) = zztmp * zce_oce(ji,jj) END DO END DO ELSE !== BLK formulation ==! turbulent fluxes computation ! ! Compute true air density : IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN ! At zu: (probably useless to remove zrho*grav*rn_zu from SLP...) zrhoa(:,:) = rho_air( t_zu(:,:) , q_zu(:,:) , pslp(:,:) ) ELSE ! At zt: zrhoa(:,:) = rho_air( ptair(:,:), phumi(:,:), pslp(:,:) ) END IF DO jj = 1, jpj DO ji = 1, jpi zztmp = zrhoa(ji,jj) * zU_zu(ji,jj) !!gm to be done by someone: check the efficiency of the call of cp_air within do loops psen (ji,jj) = cp_air( q_zu(ji,jj) ) * zztmp * zch_oce(ji,jj) * ( zst(ji,jj) - t_zu(ji,jj) ) pevp (ji,jj) = rn_efac*MAX( 0._wp, zztmp * zce_oce(ji,jj) * ( pssq(ji,jj) - q_zu(ji,jj) ) ) zztmp = zztmp * zcd_oce(ji,jj) pcd_du(ji,jj) = zztmp taum (ji,jj) = zztmp * wndm (ji,jj) zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) END DO END DO CALL iom_put( "taum_oce", taum ) ! output wind stress module ! ... utau, vtau at U- and V_points, resp. ! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines ! Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves DO jj = 1, jpjm1 DO ji = 1, fs_jpim1 utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj ) ) & & * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1)) vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji ,jj+1) ) & & * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1)) END DO END DO CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. ) IF(ln_ctl) THEN CALL prt_ctl( tab2d_1=wndm , clinfo1=' blk_oce_1: wndm : ') CALL prt_ctl( tab2d_1=utau , clinfo1=' blk_oce_1: utau : ', mask1=umask, & & tab2d_2=vtau , clinfo2=' vtau : ', mask2=vmask ) ENDIF ! ENDIF ! IF(ln_ctl) THEN CALL prt_ctl( tab2d_1=pevp , clinfo1=' blk_oce_1: pevp : ' ) CALL prt_ctl( tab2d_1=psen , clinfo1=' blk_oce_1: psen : ' ) CALL prt_ctl( tab2d_1=pssq , clinfo1=' blk_oce_1: pssq : ' ) ENDIF ! END SUBROUTINE blk_oce_1 SUBROUTINE blk_oce_2( ptair, pqsr, pqlw, pprec, & ! <<= in & psnow, pst , psen, pevp ) ! <<= in !!--------------------------------------------------------------------- !! *** ROUTINE blk_oce_2 *** !! !! ** Purpose : finalize the momentum, heat and freshwater fluxes computation !! at the ocean surface at each time step knowing Cd, Ch, Ce and !! atmospheric variables (from ABL or external data) !! !! ** Outputs : - utau : i-component of the stress at U-point (N/m2) !! - vtau : j-component of the stress at V-point (N/m2) !! - taum : Wind stress module at T-point (N/m2) !! - wndm : Wind speed module at T-point (m/s) !! - qsr : Solar heat flux over the ocean (W/m2) !! - qns : Non Solar heat flux over the ocean (W/m2) !! - emp : evaporation minus precipitation (kg/m2/s) !!--------------------------------------------------------------------- REAL(wp), INTENT(in), DIMENSION(:,:) :: ptair REAL(wp), INTENT(in), DIMENSION(:,:) :: pqsr REAL(wp), INTENT(in), DIMENSION(:,:) :: pqlw REAL(wp), INTENT(in), DIMENSION(:,:) :: pprec REAL(wp), INTENT(in), DIMENSION(:,:) :: psnow REAL(wp), INTENT(in), DIMENSION(:,:) :: pst ! surface temperature [Celcius] REAL(wp), INTENT(in), DIMENSION(:,:) :: psen REAL(wp), INTENT(in), DIMENSION(:,:) :: pevp ! INTEGER :: ji, jj ! dummy loop indices REAL(wp) :: zztmp,zz1,zz2,zz3 ! local variable REAL(wp), DIMENSION(jpi,jpj) :: zqlw ! long wave and sensible heat fluxes REAL(wp), DIMENSION(jpi,jpj) :: zqla ! latent heat fluxes and evaporation REAL(wp), DIMENSION(jpi,jpj) :: zst ! surface temperature in Kelvin !!--------------------------------------------------------------------- ! ! local scalars ( place there for vector optimisation purposes) zst(:,:) = pst(:,:) + rt0 ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K) ! ----------------------------------------------------------------------------- ! ! I Radiative FLUXES ! ! ----------------------------------------------------------------------------- ! ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle ! Short Wave zztmp = 1._wp - albo IF( ln_dm2dc ) THEN ; qsr(:,:) = zztmp * sbc_dcy( pqsr(:,:) ) * tmask(:,:,1) ELSE ; qsr(:,:) = zztmp * pqsr(:,:) * tmask(:,:,1) ENDIF zqlw(:,:) = ( pqlw(:,:) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:) ) * tmask(:,:,1) ! Long Wave ! Turbulent fluxes over ocean ! ----------------------------- zqla(:,:) = L_vap( zst(:,:) ) * pevp(:,:) ! Latent Heat flux IF(ln_ctl) THEN CALL prt_ctl( tab2d_1=zqla , clinfo1=' blk_oce_2: zqla : ' ) CALL prt_ctl( tab2d_1=zqlw , clinfo1=' blk_oce_2: zqlw : ', tab2d_2=qsr, clinfo2=' qsr : ' ) ENDIF ! ----------------------------------------------------------------------------- ! ! III Total FLUXES ! ! ----------------------------------------------------------------------------- ! ! emp (:,:) = ( pevp(:,:) - pprec(:,:) * rn_pfac ) * tmask(:,:,1) ! mass flux (evap. - precip.) ! zz1 = rn_pfac * rLfus zz2 = rn_pfac * rcp zz3 = rn_pfac * rcpi zztmp = 1._wp - albo qns(:,:) = ( zqlw(:,:) - psen(:,:) - zqla(:,:) & ! Downward Non Solar & - psnow(:,:) * zztmp & ! remove latent melting heat for solid precip & - pevp(:,:) * pst(:,:) * rcp & ! remove evap heat content at SST & + ( pprec(:,:) - psnow(:,:) ) * ( ptair(:,:) - rt0 ) * zz2 & ! add liquid precip heat content at Tair & + psnow(:,:) * ( MIN( ptair(:,:), rt0 ) - rt0 ) * zz3 & ! add solid precip heat content at min(Tair,Tsnow) & ) * tmask(:,:,1) ! #if defined key_si3 qns_oce(:,:) = zqlw(:,:) - psen(:,:) - zqla(:,:) ! non solar without emp (only needed by SI3) qsr_oce(:,:) = qsr(:,:) #endif ! IF ( nn_ice == 0 ) THEN CALL iom_put( "qlw_oce" , zqlw ) ! output downward longwave heat over the ocean CALL iom_put( "qsb_oce" , - psen ) ! output downward sensible heat over the ocean CALL iom_put( "qla_oce" , - zqla ) ! output downward latent heat over the ocean CALL iom_put( "qemp_oce", qns-zqlw+psen+zqla ) ! output downward heat content of E-P over the ocean CALL iom_put( "qns_oce" , qns ) ! output downward non solar heat over the ocean CALL iom_put( "qsr_oce" , qsr ) ! output downward solar heat over the ocean CALL iom_put( "qt_oce" , qns+qsr ) ! output total downward heat over the ocean tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s] sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s] CALL iom_put( 'snowpre', sprecip ) ! Snow CALL iom_put( 'precip' , tprecip ) ! Total precipitation ENDIF ! IF(ln_ctl) THEN CALL prt_ctl(tab2d_1=zqlw , clinfo1=' blk_oce_2: zqlw : ') CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce_2: zqla : ', tab2d_2=qsr , clinfo2=' qsr : ') CALL prt_ctl(tab2d_1=emp , clinfo1=' blk_oce_2: emp : ') ENDIF ! END SUBROUTINE blk_oce_2 FUNCTION rho_air( ptak, pqa, pslp ) !!------------------------------------------------------------------------------- !! *** FUNCTION rho_air *** !! !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere !! !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) !!------------------------------------------------------------------------------- REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K] REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg] REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp ! pressure in [Pa] REAL(wp), DIMENSION(jpi,jpj) :: rho_air ! density of moist air [kg/m^3] !!------------------------------------------------------------------------------- ! rho_air = pslp / ( R_dry*ptak * ( 1._wp + rctv0*pqa ) ) ! END FUNCTION rho_air FUNCTION cp_air_0d( pqa ) !!------------------------------------------------------------------------------- !! *** FUNCTION cp_air *** !! !! ** Purpose : Compute specific heat (Cp) of moist air !! !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) !!------------------------------------------------------------------------------- REAL(wp), INTENT(in) :: pqa ! air specific humidity [kg/kg] REAL(wp) :: cp_air_0d! specific heat of moist air [J/K/kg] !!------------------------------------------------------------------------------- ! Cp_air_0d = Cp_dry + Cp_vap * pqa ! END FUNCTION cp_air_0d FUNCTION cp_air_2d( pqa ) !!------------------------------------------------------------------------------- !! *** FUNCTION cp_air *** !! !! ** Purpose : Compute specific heat (Cp) of moist air !! !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) !!------------------------------------------------------------------------------- REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg] REAL(wp), DIMENSION(jpi,jpj) :: cp_air_2d! specific heat of moist air [J/K/kg] !!------------------------------------------------------------------------------- ! Cp_air_2d = Cp_dry + Cp_vap * pqa ! END FUNCTION cp_air_2d FUNCTION q_sat( ptak, pslp ) !!---------------------------------------------------------------------------------- !! *** FUNCTION q_sat *** !! !! ** Purpose : Specific humidity at saturation in [kg/kg] !! Based on accurate estimate of "e_sat" !! aka saturation water vapor (Goff, 1957) !! !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) !!---------------------------------------------------------------------------------- REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K] REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp ! sea level atmospheric pressure [Pa] REAL(wp), DIMENSION(jpi,jpj) :: q_sat ! Specific humidity at saturation [kg/kg] ! INTEGER :: ji, jj ! dummy loop indices REAL(wp) :: ze_sat, ztmp ! local scalar !!---------------------------------------------------------------------------------- ! DO jj = 1, jpj DO ji = 1, jpi ! ztmp = rt0 / ptak(ji,jj) ! ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957) ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0) & & + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) ) & & + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614 ) ! q_sat(ji,jj) = reps0 * ze_sat/( 0.01_wp*pslp(ji,jj) - (1._wp - reps0)*ze_sat ) ! 0.01 because SLP is in [Pa] ! END DO END DO ! END FUNCTION q_sat FUNCTION gamma_moist( ptak, pqa ) !!---------------------------------------------------------------------------------- !! *** FUNCTION gamma_moist *** !! !! ** Purpose : Compute the moist adiabatic lapse-rate. !! => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate !! => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html !! !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) !!---------------------------------------------------------------------------------- REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K] REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! specific humidity [kg/kg] REAL(wp), DIMENSION(jpi,jpj) :: gamma_moist ! moist adiabatic lapse-rate ! INTEGER :: ji, jj ! dummy loop indices REAL(wp) :: zrv, ziRT ! local scalar !!---------------------------------------------------------------------------------- ! DO jj = 1, jpj DO ji = 1, jpi zrv = pqa(ji,jj) / (1. - pqa(ji,jj)) ziRT = 1. / (R_dry*ptak(ji,jj)) ! 1/RT gamma_moist(ji,jj) = grav * ( 1. + rLevap*zrv*ziRT ) / ( Cp_dry + rLevap*rLevap*zrv*reps0*ziRT/ptak(ji,jj) ) END DO END DO ! END FUNCTION gamma_moist FUNCTION L_vap( psst ) !!--------------------------------------------------------------------------------- !! *** FUNCTION L_vap *** !! !! ** Purpose : Compute the latent heat of vaporization of water from temperature !! !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) !!---------------------------------------------------------------------------------- REAL(wp), DIMENSION(jpi,jpj) :: L_vap ! latent heat of vaporization [J/kg] REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst ! water temperature [K] !!---------------------------------------------------------------------------------- ! L_vap = ( 2.501 - 0.00237 * ( psst(:,:) - rt0) ) * 1.e6 ! END FUNCTION L_vap #if defined key_si3 !!---------------------------------------------------------------------- !! 'key_si3' SI3 sea-ice model !!---------------------------------------------------------------------- !! blk_ice_1 : provide the air-ice stress !! blk_ice_2 : provide the heat and mass fluxes at air-ice interface !! blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux) !! Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag !! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag !!---------------------------------------------------------------------- SUBROUTINE blk_ice_1( pwndi, pwndj, ptair, phumi, pslp , puice, pvice, ptsui, & ! inputs & putaui, pvtaui, pseni, pevpi, pssqi, pcd_dui ) ! optional outputs !!--------------------------------------------------------------------- !! *** ROUTINE blk_ice_1 *** !! !! ** Purpose : provide the surface boundary condition over sea-ice !! !! ** Method : compute momentum using bulk formulation !! formulea, ice variables and read atmospheric fields. !! NB: ice drag coefficient is assumed to be a constant !!--------------------------------------------------------------------- REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pslp ! sea-level pressure [Pa] REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndi ! atmospheric wind at T-point [m/s] REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndj ! atmospheric wind at T-point [m/s] REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: ptair ! atmospheric wind at T-point [m/s] REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: phumi ! atmospheric wind at T-point [m/s] REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: puice ! sea-ice velocity on I or C grid [m/s] REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pvice ! " REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: ptsui ! sea-ice surface temperature [K] REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: putaui ! if ln_blk REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pvtaui ! if ln_blk REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pseni ! if ln_abl REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pevpi ! if ln_abl REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pssqi ! if ln_abl REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pcd_dui ! if ln_abl ! INTEGER :: ji, jj ! dummy loop indices REAL(wp) :: zwndi_t , zwndj_t ! relative wind components at T-point REAL(wp) :: zootm_su ! sea-ice surface mean temperature REAL(wp) :: zztmp1, zztmp2 ! temporary arrays REAL(wp), DIMENSION(jpi,jpj) :: zrhoa ! transfer coefficient for momentum (tau) REAL(wp), DIMENSION(jpi,jpj) :: zcd_dui ! transfer coefficient for momentum (tau) !!--------------------------------------------------------------------- ! ------------------------------------------------------------ ! ! Wind module relative to the moving ice ( U10m - U_ice ) ! ! ------------------------------------------------------------ ! ! C-grid ice dynamics : U & V-points (same as ocean) DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vect. opt. zwndi_t = ( pwndi(ji,jj) - rn_vfac * 0.5_wp * ( puice(ji-1,jj ) + puice(ji,jj) ) ) zwndj_t = ( pwndj(ji,jj) - rn_vfac * 0.5_wp * ( pvice(ji ,jj-1) + pvice(ji,jj) ) ) wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) END DO END DO CALL lbc_lnk( 'sbcblk', wndm_ice, 'T', 1. ) ! ! Make ice-atm. drag dependent on ice concentration IF ( ln_Cd_L12 ) THEN ! calculate new drag from Lupkes(2012) equations CALL Cdn10_Lupkes2012( Cd_ice ) Ch_ice(:,:) = Cd_ice(:,:) ! momentum and heat transfer coef. are considered identical Ce_ice(:,:) = Cd_ice(:,:) ELSEIF( ln_Cd_L15 ) THEN ! calculate new drag from Lupkes(2015) equations CALL Cdn10_Lupkes2015( ptsui, Cd_ice, Ch_ice ) Ce_ice(:,:) = Ch_ice(:,:) ! sensible and latent heat transfer coef. are considered identical ENDIF IF ( iom_use("Cd_ice") ) CALL iom_put("Cd_ice", Cd_ice) ! output value of pure ice-atm. transfer coef. IF ( iom_use("Ch_ice") ) CALL iom_put("Ch_ice", Ch_ice) ! output value of pure ice-atm. transfer coef. ! local scalars ( place there for vector optimisation purposes) ! Computing density of air! Way denser that 1.2 over sea-ice !!! zrhoa (:,:) = rho_air( ptair(:,:), phumi(:,:), pslp(:,:) ) zcd_dui(:,:) = wndm_ice(:,:) * Cd_ice(:,:) IF( ln_blk ) THEN ! ------------------------------------------------------------ ! ! Wind stress relative to the moving ice ( U10m - U_ice ) ! ! ------------------------------------------------------------ ! ! C-grid ice dynamics : U & V-points (same as ocean) DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vect. opt. putaui(ji,jj) = 0.5_wp * ( zrhoa(ji+1,jj) * zcd_dui(ji+1,jj) & & + zrhoa(ji ,jj) * zcd_dui(ji ,jj) ) & & * ( 0.5_wp * ( pwndi(ji+1,jj) + pwndi(ji,jj) ) - rn_vfac * puice(ji,jj) ) pvtaui(ji,jj) = 0.5_wp * ( zrhoa(ji,jj+1) * zcd_dui(ji,jj+1) & & + zrhoa(ji,jj ) * zcd_dui(ji,jj ) ) & & * ( 0.5_wp * ( pwndj(ji,jj+1) + pwndj(ji,jj) ) - rn_vfac * pvice(ji,jj) ) END DO END DO CALL lbc_lnk_multi( 'sbcblk', putaui, 'U', -1., pvtaui, 'V', -1. ) ! IF(ln_ctl) CALL prt_ctl( tab2d_1=putaui , clinfo1=' blk_ice: putaui : ' & & , tab2d_2=pvtaui , clinfo2=' pvtaui : ' ) ELSE zztmp1 = 11637800.0_wp zztmp2 = -5897.8_wp DO jj = 1, jpj DO ji = 1, jpi pcd_dui(ji,jj) = zcd_dui (ji,jj) pseni (ji,jj) = wndm_ice(ji,jj) * Ch_ice(ji,jj) pevpi (ji,jj) = wndm_ice(ji,jj) * Ce_ice(ji,jj) zootm_su = zztmp2 / ptsui(ji,jj) ! ptsui is in K (it can't be zero ??) pssqi (ji,jj) = zztmp1 * EXP( zootm_su ) / zrhoa(ji,jj) END DO END DO ENDIF ! IF(ln_ctl) CALL prt_ctl(tab2d_1=wndm_ice , clinfo1=' blk_ice: wndm_ice : ') ! END SUBROUTINE blk_ice_1 SUBROUTINE blk_ice_2( ptsu, phs, phi, palb, ptair, phumi, pslp, pqlw, pprec, psnow ) !!--------------------------------------------------------------------- !! *** ROUTINE blk_ice_2 *** !! !! ** Purpose : provide the heat and mass fluxes at air-ice interface !! !! ** Method : compute heat and freshwater exchanged !! between atmosphere and sea-ice using bulk formulation !! formulea, ice variables and read atmmospheric fields. !! !! caution : the net upward water flux has with mm/day unit !!--------------------------------------------------------------------- REAL(wp), DIMENSION(:,:,:), INTENT(in) :: ptsu ! sea ice surface temperature [K] REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phs ! snow thickness REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phi ! ice thickness REAL(wp), DIMENSION(:,:,:), INTENT(in) :: palb ! ice albedo (all skies) REAL(wp), DIMENSION(:,: ), INTENT(in) :: ptair REAL(wp), DIMENSION(:,: ), INTENT(in) :: phumi REAL(wp), DIMENSION(:,: ), INTENT(in) :: pslp REAL(wp), DIMENSION(:,: ), INTENT(in) :: pqlw REAL(wp), DIMENSION(:,: ), INTENT(in) :: pprec REAL(wp), DIMENSION(:,: ), INTENT(in) :: psnow !! INTEGER :: ji, jj, jl ! dummy loop indices REAL(wp) :: zst3 ! local variable REAL(wp) :: zcoef_dqlw, zcoef_dqla ! - - REAL(wp) :: zztmp, zztmp2, z1_rLsub ! - - REAL(wp) :: zfr1, zfr2 ! local variables REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_st ! inverse of surface temperature REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qlw ! long wave heat flux over ice REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_qsb ! sensible heat flux over ice REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqlw ! long wave heat sensitivity over ice REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqsb ! sensible heat sensitivity over ice REAL(wp), DIMENSION(jpi,jpj) :: zevap, zsnw ! evaporation and snw distribution after wind blowing (SI3) REAL(wp), DIMENSION(jpi,jpj) :: zrhoa !!--------------------------------------------------------------------- ! zcoef_dqlw = 4.0 * 0.95 * Stef ! local scalars zcoef_dqla = -Ls * 11637800. * (-5897.8) ! zrhoa(:,:) = rho_air( ptair(:,:), phumi(:,:), pslp(:,:) ) ! zztmp = 1. / ( 1. - albo ) WHERE( ptsu(:,:,:) /= 0._wp ) ; z1_st(:,:,:) = 1._wp / ptsu(:,:,:) ELSEWHERE ; z1_st(:,:,:) = 0._wp END WHERE ! ! ========================== ! DO jl = 1, jpl ! Loop over ice categories ! ! ! ========================== ! DO jj = 1 , jpj DO ji = 1, jpi ! ----------------------------! ! I Radiative FLUXES ! ! ----------------------------! zst3 = ptsu(ji,jj,jl) * ptsu(ji,jj,jl) * ptsu(ji,jj,jl) ! Short Wave (sw) qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj) ! Long Wave (lw) z_qlw(ji,jj,jl) = 0.95 * ( pqlw(ji,jj) - Stef * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) ! lw sensitivity z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3 ! ----------------------------! ! II Turbulent FLUXES ! ! ----------------------------! ! ... turbulent heat fluxes with Ch_ice recalculated in blk_ice_1 ! Sensible Heat z_qsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_ice(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - ptair(ji,jj)) ! Latent Heat zztmp2 = EXP( -5897.8 * z1_st(ji,jj,jl) ) qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls * Ce_ice(ji,jj) * wndm_ice(ji,jj) * & & ( 11637800. * zztmp2 / zrhoa(ji,jj) - phumi(ji,jj) ) ) ! Latent heat sensitivity for ice (Dqla/Dt) IF( qla_ice(ji,jj,jl) > 0._wp ) THEN dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ce_ice(ji,jj) * wndm_ice(ji,jj) * & & z1_st(ji,jj,jl) * z1_st(ji,jj,jl) * zztmp2 ELSE dqla_ice(ji,jj,jl) = 0._wp ENDIF ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice) z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_ice(ji,jj) * wndm_ice(ji,jj) ! ----------------------------! ! III Total FLUXES ! ! ----------------------------! ! Downward Non Solar flux qns_ice (ji,jj,jl) = z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl) ! Total non solar heat flux sensitivity for ice dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) ) END DO ! END DO ! END DO ! tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1) ! total precipitation [kg/m2/s] sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1) ! solid precipitation [kg/m2/s] CALL iom_put( 'snowpre', sprecip ) ! Snow precipitation CALL iom_put( 'precip' , tprecip ) ! Total precipitation ! --- evaporation --- ! z1_rLsub = 1._wp / rLsub evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_rLsub ! sublimation devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_rLsub ! d(sublimation)/dT zevap (:,:) = rn_efac * ( emp(:,:) + tprecip(:,:) ) ! evaporation over ocean ! --- evaporation minus precipitation --- ! zsnw(:,:) = 0._wp CALL ice_thd_snwblow( (1.-at_i_b(:,:)), zsnw ) ! snow distribution over ice after wind blowing emp_oce(:,:) = ( 1._wp - at_i_b(:,:) ) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw ) emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:) ! --- heat flux associated with emp --- ! qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp & ! evap at sst & + ( tprecip(:,:) - sprecip(:,:) ) * ( ptair(:,:) - rt0 ) * rcp & ! liquid precip at Tair & + sprecip(:,:) * ( 1._wp - zsnw ) * & ! solid precip at min(Tair,Tsnow) & ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) qemp_ice(:,:) = sprecip(:,:) * zsnw * & ! solid precip (only) & ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) ! --- total solar and non solar fluxes --- ! qns_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) & & + qemp_ice(:,:) + qemp_oce(:,:) qsr_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 ) ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- ! qprec_ice(:,:) = rhos * ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- DO jl = 1, jpl qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * rcpi * tmask(:,:,1) ) ! ! But we do not have Tice => consider it at 0degC => evap=0 END DO ! --- shortwave radiation transmitted below the surface (W/m2, see Grenfell Maykut 77) --- ! zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) ! transmission when hi>10cm zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) ! zfr2 such that zfr1 + zfr2 to equal 1 ! WHERE ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) < 0.1_wp ) ! linear decrease from hi=0 to 10cm qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) ) ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp ) ! constant (zfr1) when hi>10cm qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1 ELSEWHERE ! zero when hs>0 qtr_ice_top(:,:,:) = 0._wp END WHERE ! IF(ln_ctl) THEN CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice : ', tab3d_2=z_qsb , clinfo2=' z_qsb : ', kdim=jpl) CALL prt_ctl(tab3d_1=z_qlw , clinfo1=' blk_ice: z_qlw : ', tab3d_2=dqla_ice, clinfo2=' dqla_ice : ', kdim=jpl) CALL prt_ctl(tab3d_1=z_dqsb , clinfo1=' blk_ice: z_dqsb : ', tab3d_2=z_dqlw , clinfo2=' z_dqlw : ', kdim=jpl) CALL prt_ctl(tab3d_1=dqns_ice, clinfo1=' blk_ice: dqns_ice : ', tab3d_2=qsr_ice , clinfo2=' qsr_ice : ', kdim=jpl) CALL prt_ctl(tab3d_1=ptsu , clinfo1=' blk_ice: ptsu : ', tab3d_2=qns_ice , clinfo2=' qns_ice : ', kdim=jpl) CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice: tprecip : ', tab2d_2=sprecip , clinfo2=' sprecip : ') END IF ! END SUBROUTINE blk_ice_2 SUBROUTINE blk_ice_qcn( ld_virtual_itd, ptsu, ptb, phs, phi ) !!--------------------------------------------------------------------- !! *** ROUTINE blk_ice_qcn *** !! !! ** Purpose : Compute surface temperature and snow/ice conduction flux !! to force sea ice / snow thermodynamics !! in the case conduction flux is emulated !! !! ** Method : compute surface energy balance assuming neglecting heat storage !! following the 0-layer Semtner (1976) approach !! !! ** Outputs : - ptsu : sea-ice / snow surface temperature (K) !! - qcn_ice : surface inner conduction flux (W/m2) !! !!--------------------------------------------------------------------- LOGICAL , INTENT(in ) :: ld_virtual_itd ! single-category option REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: ptsu ! sea ice / snow surface temperature REAL(wp), DIMENSION(:,:) , INTENT(in ) :: ptb ! sea ice base temperature REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: phs ! snow thickness REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: phi ! sea ice thickness ! INTEGER , PARAMETER :: nit = 10 ! number of iterations REAL(wp), PARAMETER :: zepsilon = 0.1_wp ! characteristic thickness for enhanced conduction ! INTEGER :: ji, jj, jl ! dummy loop indices INTEGER :: iter ! local integer REAL(wp) :: zfac, zfac2, zfac3 ! local scalars REAL(wp) :: zkeff_h, ztsu, ztsu0 ! REAL(wp) :: zqc, zqnet ! REAL(wp) :: zhe, zqa0 ! REAL(wp), DIMENSION(jpi,jpj,jpl) :: zgfac ! enhanced conduction factor !!--------------------------------------------------------------------- ! -------------------------------------! ! I Enhanced conduction factor ! ! -------------------------------------! ! Emulates the enhancement of conduction by unresolved thin ice (ld_virtual_itd = T) ! Fichefet and Morales Maqueda, JGR 1997 ! zgfac(:,:,:) = 1._wp IF( ld_virtual_itd ) THEN ! zfac = 1._wp / ( rn_cnd_s + rcnd_i ) zfac2 = EXP(1._wp) * 0.5_wp * zepsilon zfac3 = 2._wp / zepsilon ! DO jl = 1, jpl DO jj = 1 , jpj DO ji = 1, jpi zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac ! Effective thickness IF( zhe >= zfac2 ) zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor END DO END DO END DO ! ENDIF ! -------------------------------------------------------------! ! II Surface temperature and conduction flux ! ! -------------------------------------------------------------! ! zfac = rcnd_i * rn_cnd_s ! DO jl = 1, jpl DO jj = 1 , jpj DO ji = 1, jpi ! zkeff_h = zfac * zgfac(ji,jj,jl) / & ! Effective conductivity of the snow-ice system divided by thickness & ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) ) ztsu = ptsu(ji,jj,jl) ! Store current iteration temperature ztsu0 = ptsu(ji,jj,jl) ! Store initial surface temperature zqa0 = qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) ! Net initial atmospheric heat flux ! DO iter = 1, nit ! --- Iterative loop zqc = zkeff_h * ( ztsu - ptb(ji,jj) ) ! Conduction heat flux through snow-ice system (>0 downwards) zqnet = zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc ! Surface energy budget ztsu = ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h ) ! Temperature update END DO ! ptsu (ji,jj,jl) = MIN( rt0, ztsu ) qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) ) qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) qml_ice(ji,jj,jl) = ( qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) - qcn_ice(ji,jj,jl) ) & & * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) ) ! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- ! hfx_err_dif(ji,jj) = hfx_err_dif(ji,jj) - ( dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) ) * a_i_b(ji,jj,jl) END DO END DO ! END DO ! END SUBROUTINE blk_ice_qcn SUBROUTINE Cdn10_Lupkes2012( pcd ) !!---------------------------------------------------------------------- !! *** ROUTINE Cdn10_Lupkes2012 *** !! !! ** Purpose : Recompute the neutral air-ice drag referenced at 10m !! to make it dependent on edges at leads, melt ponds and flows. !! After some approximations, this can be resumed to a dependency !! on ice concentration. !! !! ** Method : The parameterization is taken from Lupkes et al. (2012) eq.(50) !! with the highest level of approximation: level4, eq.(59) !! The generic drag over a cell partly covered by ice can be re-written as follows: !! !! Cd = Cdw * (1-A) + Cdi * A + Ce * (1-A)**(nu+1/(10*beta)) * A**mu !! !! Ce = 2.23e-3 , as suggested by Lupkes (eq. 59) !! nu = mu = beta = 1 , as suggested by Lupkes (eq. 59) !! A is the concentration of ice minus melt ponds (if any) !! !! This new drag has a parabolic shape (as a function of A) starting at !! Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5 !! and going down to Cdi(say 1.4e-3) for A=1 !! !! It is theoretically applicable to all ice conditions (not only MIZ) !! => see Lupkes et al (2013) !! !! ** References : Lupkes et al. JGR 2012 (theory) !! Lupkes et al. GRL 2013 (application to GCM) !! !!---------------------------------------------------------------------- REAL(wp), DIMENSION(:,:), INTENT(inout) :: pcd REAL(wp), PARAMETER :: zCe = 2.23e-03_wp REAL(wp), PARAMETER :: znu = 1._wp REAL(wp), PARAMETER :: zmu = 1._wp REAL(wp), PARAMETER :: zbeta = 1._wp REAL(wp) :: zcoef !!---------------------------------------------------------------------- zcoef = znu + 1._wp / ( 10._wp * zbeta ) ! generic drag over a cell partly covered by ice !!Cd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) + & ! pure ocean drag !! & Cd_ice * at_i_b(:,:) + & ! pure ice drag !! & zCe * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**zmu ! change due to sea-ice morphology ! ice-atm drag pcd(:,:) = rcdice + & ! pure ice drag & zCe * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp) ! change due to sea-ice morphology END SUBROUTINE Cdn10_Lupkes2012 SUBROUTINE Cdn10_Lupkes2015( ptm_su, pcd, pch ) !!---------------------------------------------------------------------- !! *** ROUTINE Cdn10_Lupkes2015 *** !! !! ** pUrpose : Alternative turbulent transfert coefficients formulation !! between sea-ice and atmosphere with distinct momentum !! and heat coefficients depending on sea-ice concentration !! and atmospheric stability (no meltponds effect for now). !! !! ** Method : The parameterization is adapted from Lupkes et al. (2015) !! and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme, !! it considers specific skin and form drags (Andreas et al. 2010) !! to compute neutral transfert coefficients for both heat and !! momemtum fluxes. Atmospheric stability effect on transfert !! coefficient is also taken into account following Louis (1979). !! !! ** References : Lupkes et al. JGR 2015 (theory) !! Lupkes et al. ECHAM6 documentation 2015 (implementation) !! !!---------------------------------------------------------------------- ! REAL(wp), DIMENSION(:,:), INTENT(in ) :: ptm_su ! sea-ice surface temperature [K] REAL(wp), DIMENSION(:,:), INTENT(inout) :: pcd ! momentum transfert coefficient REAL(wp), DIMENSION(:,:), INTENT(inout) :: pch ! heat transfert coefficient REAL(wp), DIMENSION(jpi,jpj) :: zst, zqo_sat, zqi_sat ! ! ECHAM6 constants REAL(wp), PARAMETER :: z0_skin_ice = 0.69e-3_wp ! Eq. 43 [m] REAL(wp), PARAMETER :: z0_form_ice = 0.57e-3_wp ! Eq. 42 [m] REAL(wp), PARAMETER :: z0_ice = 1.00e-3_wp ! Eq. 15 [m] REAL(wp), PARAMETER :: zce10 = 2.80e-3_wp ! Eq. 41 REAL(wp), PARAMETER :: zbeta = 1.1_wp ! Eq. 41 REAL(wp), PARAMETER :: zc = 5._wp ! Eq. 13 REAL(wp), PARAMETER :: zc2 = zc * zc REAL(wp), PARAMETER :: zam = 2. * zc ! Eq. 14 REAL(wp), PARAMETER :: zah = 3. * zc ! Eq. 30 REAL(wp), PARAMETER :: z1_alpha = 1._wp / 0.2_wp ! Eq. 51 REAL(wp), PARAMETER :: z1_alphaf = z1_alpha ! Eq. 56 REAL(wp), PARAMETER :: zbetah = 1.e-3_wp ! Eq. 26 REAL(wp), PARAMETER :: zgamma = 1.25_wp ! Eq. 26 REAL(wp), PARAMETER :: z1_gamma = 1._wp / zgamma REAL(wp), PARAMETER :: r1_3 = 1._wp / 3._wp ! INTEGER :: ji, jj ! dummy loop indices REAL(wp) :: zthetav_os, zthetav_is, zthetav_zu REAL(wp) :: zrib_o, zrib_i REAL(wp) :: zCdn_skin_ice, zCdn_form_ice, zCdn_ice REAL(wp) :: zChn_skin_ice, zChn_form_ice REAL(wp) :: z0w, z0i, zfmi, zfmw, zfhi, zfhw REAL(wp) :: zCdn_form_tmp !!--------------------------------------------------------------------------- ! ! Momentum Neutral Transfert Coefficients (should be a constant) zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2 ! Eq. 40 zCdn_skin_ice = ( vkarmn / LOG( rn_zu / z0_skin_ice + 1._wp ) )**2 ! Eq. 7 zCdn_ice = zCdn_skin_ice ! Eq. 7 !zCdn_ice = 1.89e-3 ! old ECHAM5 value (cf Eq. 32) ! Heat Neutral Transfert Coefficients 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 ! Atmospheric and Surface Variables zst(:,:) = sst_m(:,:) + rt0 ! convert SST from Celcius to Kelvin zqo_sat(:,:) = 0.98_wp * q_sat( zst(:,:) , sf(jp_slp)%fnow(:,:,1) ) ! saturation humidity over ocean [kg/kg] zqi_sat(:,:) = 0.98_wp * q_sat( ptm_su(:,:), sf(jp_slp)%fnow(:,:,1) ) ! saturation humidity over ice [kg/kg] ! DO jj = 2, jpjm1 ! reduced loop is necessary for reproductibility DO ji = fs_2, fs_jpim1 ! Virtual potential temperature [K] zthetav_os = zst(ji,jj) * ( 1._wp + rctv0 * zqo_sat(ji,jj) ) ! over ocean zthetav_is = ptm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) ) ! ocean ice zthetav_zu = t_zu (ji,jj) * ( 1._wp + rctv0 * q_zu(ji,jj) ) ! at zu ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead) zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj) )**2 ! over ocean zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2 ! over ice ! Momentum and Heat Neutral Transfert Coefficients zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta ! Eq. 40 zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) ) ! Eq. 53 ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead) z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water z0i = z0_skin_ice ! over ice IF( zrib_o <= 0._wp ) THEN 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 zfhw = ( 1._wp + ( zbetah * ( zthetav_os - zthetav_zu )**r1_3 / ( Chn_oce(ji,jj) * MAX(0.01, wndm(ji,jj)) ) & ! Eq. 26 & )**zgamma )**z1_gamma ELSE zfmw = 1._wp / ( 1._wp + zam * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 12 zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 28 ENDIF IF( zrib_i <= 0._wp ) THEN zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 9 zfhi = 1._wp - zah * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 25 ELSE zfmi = 1._wp / ( 1._wp + zam * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 11 zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 27 ENDIF ! Momentum Transfert Coefficients (Eq. 38) pcd(ji,jj) = zCdn_skin_ice * zfmi + & & 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) ) ! Heat Transfert Coefficients (Eq. 49) pch(ji,jj) = zChn_skin_ice * zfhi + & & 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) ) ! END DO END DO ! CALL lbc_lnk_multi( 'sbcblk', pcd, 'T', 1., pch, 'T', 1. ) ! END SUBROUTINE Cdn10_Lupkes2015 #endif !!====================================================================== END MODULE sbcblk