- Timestamp:
- 2020-09-14T17:40:34+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/SBC/sbcblk.F90 (modified) (40 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEADext/AGRIF5 ^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@13382 sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/OCE/SBC/sbcblk.F90
r11482 r13463 15 15 !! 3.7 ! 2014-06 (L. Brodeau) simplification and optimization of CORE bulk 16 16 !! 4.0 ! 2016-06 (L. Brodeau) sbcblk_core becomes sbcblk and is not restricted to the CORE algorithm anymore 17 !! ! ==> based on AeroBulk (http ://aerobulk.sourceforge.net/)17 !! ! ==> based on AeroBulk (https://github.com/brodeau/aerobulk/) 18 18 !! 4.0 ! 2016-10 (G. Madec) introduce a sbc_blk_init routine 19 !! 4.0 ! 2016-10 (M. Vancoppenolle) Introduce conduction flux emulator (M. Vancoppenolle) 19 !! 4.0 ! 2016-10 (M. Vancoppenolle) Introduce conduction flux emulator (M. Vancoppenolle) 20 !! 4.0 ! 2019-03 (F. Lemarié & G. Samson) add ABL compatibility (ln_abl=TRUE) 20 21 !!---------------------------------------------------------------------- 21 22 … … 23 24 !! sbc_blk_init : initialisation of the chosen bulk formulation as ocean surface boundary condition 24 25 !! sbc_blk : bulk formulation as ocean surface boundary condition 25 !! blk_oce : computes momentum, heat and freshwater fluxes over ocean 26 !! rho_air : density of (moist) air (depends on T_air, q_air and SLP 27 !! cp_air : specific heat of (moist) air (depends spec. hum. q_air) 28 !! q_sat : saturation humidity as a function of SLP and temperature 29 !! L_vap : latent heat of vaporization of water as a function of temperature 30 !! sea-ice case only : 31 !! blk_ice_tau : provide the air-ice stress 32 !! blk_ice_flx : provide the heat and mass fluxes at air-ice interface 26 !! blk_oce_1 : computes pieces of momentum, heat and freshwater fluxes over ocean for ABL model (ln_abl=TRUE) 27 !! blk_oce_2 : finalizes momentum, heat and freshwater fluxes computation over ocean after the ABL step (ln_abl=TRUE) 28 !! sea-ice case only : 29 !! blk_ice_1 : provide the air-ice stress 30 !! blk_ice_2 : provide the heat and mass fluxes at air-ice interface 33 31 !! blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux) 34 32 !! Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag 35 !! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 33 !! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 36 34 !!---------------------------------------------------------------------- 37 35 USE oce ! ocean dynamics and tracers … … 46 44 USE lib_fortran ! to use key_nosignedzero 47 45 #if defined key_si3 48 USE ice , ONLY : u_ice, v_ice, jpl, a_i_b, at_i_b, t_su, rn_cnd_s, hfx_err_dif46 USE ice , ONLY : jpl, a_i_b, at_i_b, rn_cnd_s, hfx_err_dif 49 47 USE icethd_dh ! for CALL ice_thd_snwblow 50 48 #endif 51 USE sbcblk_algo_ncar ! => turb_ncar : NCAR - CORE (Large & Yeager, 2009) 52 USE sbcblk_algo_coare ! => turb_coare : COAREv3.0 (Fairall et al. 2003)53 USE sbcblk_algo_coare3p 5 ! => turb_coare3p5 : COAREv3.5 (Edson et al. 2013)54 USE sbcblk_algo_ecmwf ! => turb_ecmwf : ECMWF (IFS cycle 31)49 USE sbcblk_algo_ncar ! => turb_ncar : NCAR - CORE (Large & Yeager, 2009) 50 USE sbcblk_algo_coare3p0 ! => turb_coare3p0 : COAREv3.0 (Fairall et al. 2003) 51 USE sbcblk_algo_coare3p6 ! => turb_coare3p6 : COAREv3.6 (Fairall et al. 2018 + Edson et al. 2013) 52 USE sbcblk_algo_ecmwf ! => turb_ecmwf : ECMWF (IFS cycle 45r1) 55 53 ! 56 54 USE iom ! I/O manager library … … 60 58 USE prtctl ! Print control 61 59 60 USE sbcblk_phy ! a catalog of functions for physical/meteorological parameters in the marine boundary layer, rho_air, q_sat, etc... 61 62 62 63 IMPLICIT NONE 63 64 PRIVATE … … 65 66 PUBLIC sbc_blk_init ! called in sbcmod 66 67 PUBLIC sbc_blk ! called in sbcmod 68 PUBLIC blk_oce_1 ! called in sbcabl 69 PUBLIC blk_oce_2 ! called in sbcabl 67 70 #if defined key_si3 68 PUBLIC blk_ice_ tau! routine called in icesbc69 PUBLIC blk_ice_ flx! routine called in icesbc71 PUBLIC blk_ice_1 ! routine called in icesbc 72 PUBLIC blk_ice_2 ! routine called in icesbc 70 73 PUBLIC blk_ice_qcn ! routine called in icesbc 71 #endif 72 73 !!Lolo: should ultimately be moved in the module with all physical constants ? 74 !!gm : In principle, yes. 75 REAL(wp), PARAMETER :: Cp_dry = 1005.0 !: Specic heat of dry air, constant pressure [J/K/kg] 76 REAL(wp), PARAMETER :: Cp_vap = 1860.0 !: Specic heat of water vapor, constant pressure [J/K/kg] 77 REAL(wp), PARAMETER :: R_dry = 287.05_wp !: Specific gas constant for dry air [J/K/kg] 78 REAL(wp), PARAMETER :: R_vap = 461.495_wp !: Specific gas constant for water vapor [J/K/kg] 79 REAL(wp), PARAMETER :: reps0 = R_dry/R_vap !: ratio of gas constant for dry air and water vapor => ~ 0.622 80 REAL(wp), PARAMETER :: rctv0 = R_vap/R_dry !: for virtual temperature (== (1-eps)/eps) => ~ 0.608 81 82 INTEGER , PARAMETER :: jpfld =10 ! maximum number of files to read 83 INTEGER , PARAMETER :: jp_wndi = 1 ! index of 10m wind velocity (i-component) (m/s) at T-point 84 INTEGER , PARAMETER :: jp_wndj = 2 ! index of 10m wind velocity (j-component) (m/s) at T-point 85 INTEGER , PARAMETER :: jp_tair = 3 ! index of 10m air temperature (Kelvin) 86 INTEGER , PARAMETER :: jp_humi = 4 ! index of specific humidity ( % ) 87 INTEGER , PARAMETER :: jp_qsr = 5 ! index of solar heat (W/m2) 88 INTEGER , PARAMETER :: jp_qlw = 6 ! index of Long wave (W/m2) 89 INTEGER , PARAMETER :: jp_prec = 7 ! index of total precipitation (rain+snow) (Kg/m2/s) 90 INTEGER , PARAMETER :: jp_snow = 8 ! index of snow (solid prcipitation) (kg/m2/s) 91 INTEGER , PARAMETER :: jp_slp = 9 ! index of sea level pressure (Pa) 92 INTEGER , PARAMETER :: jp_tdif =10 ! index of tau diff associated to HF tau (N/m2) at T-point 93 94 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf ! structure of input fields (file informations, fields read) 95 96 ! !!! Bulk parameters 97 REAL(wp), PARAMETER :: cpa = 1000.5 ! specific heat of air (only used for ice fluxes now...) 98 REAL(wp), PARAMETER :: Ls = 2.839e6 ! latent heat of sublimation 99 REAL(wp), PARAMETER :: Stef = 5.67e-8 ! Stefan Boltzmann constant 100 REAL(wp), PARAMETER :: Cd_ice = 1.4e-3 ! transfer coefficient over ice 101 REAL(wp), PARAMETER :: albo = 0.066 ! ocean albedo assumed to be constant 102 ! 74 #endif 75 76 INTEGER , PUBLIC, PARAMETER :: jp_wndi = 1 ! index of 10m wind velocity (i-component) (m/s) at T-point 77 INTEGER , PUBLIC, PARAMETER :: jp_wndj = 2 ! index of 10m wind velocity (j-component) (m/s) at T-point 78 INTEGER , PUBLIC, PARAMETER :: jp_tair = 3 ! index of 10m air temperature (Kelvin) 79 INTEGER , PUBLIC, PARAMETER :: jp_humi = 4 ! index of specific humidity ( % ) 80 INTEGER , PUBLIC, PARAMETER :: jp_qsr = 5 ! index of solar heat (W/m2) 81 INTEGER , PUBLIC, PARAMETER :: jp_qlw = 6 ! index of Long wave (W/m2) 82 INTEGER , PUBLIC, PARAMETER :: jp_prec = 7 ! index of total precipitation (rain+snow) (Kg/m2/s) 83 INTEGER , PUBLIC, PARAMETER :: jp_snow = 8 ! index of snow (solid prcipitation) (kg/m2/s) 84 INTEGER , PUBLIC, PARAMETER :: jp_slp = 9 ! index of sea level pressure (Pa) 85 INTEGER , PUBLIC, PARAMETER :: jp_uoatm = 10 ! index of surface current (i-component) 86 ! ! seen by the atmospheric forcing (m/s) at T-point 87 INTEGER , PUBLIC, PARAMETER :: jp_voatm = 11 ! index of surface current (j-component) 88 ! ! seen by the atmospheric forcing (m/s) at T-point 89 INTEGER , PUBLIC, PARAMETER :: jp_hpgi = 12 ! index of ABL geostrophic wind or hpg (i-component) (m/s) at T-point 90 INTEGER , PUBLIC, PARAMETER :: jp_hpgj = 13 ! index of ABL geostrophic wind or hpg (j-component) (m/s) at T-point 91 INTEGER , PUBLIC, PARAMETER :: jpfld = 13 ! maximum number of files to read 92 93 ! Warning: keep this structure allocatable for Agrif... 94 TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) :: sf ! structure of input atmospheric fields (file informations, fields read) 95 103 96 ! !!* Namelist namsbc_blk : bulk parameters 104 97 LOGICAL :: ln_NCAR ! "NCAR" algorithm (Large and Yeager 2008) 105 98 LOGICAL :: ln_COARE_3p0 ! "COARE 3.0" algorithm (Fairall et al. 2003) 106 LOGICAL :: ln_COARE_3p 5 ! "COARE 3.5" algorithm (Edson et al. 2013)107 LOGICAL :: ln_ECMWF ! "ECMWF" algorithm (IFS cycle 31)99 LOGICAL :: ln_COARE_3p6 ! "COARE 3.6" algorithm (Edson et al. 2013) 100 LOGICAL :: ln_ECMWF ! "ECMWF" algorithm (IFS cycle 45r1) 108 101 ! 109 LOGICAL :: ln_taudif ! logical flag to use the "mean of stress module - module of mean stress" data 110 REAL(wp) :: rn_pfac ! multiplication factor for precipitation 111 REAL(wp) :: rn_efac ! multiplication factor for evaporation 112 REAL(wp) :: rn_vfac ! multiplication factor for ice/ocean velocity in the calculation of wind stress 113 REAL(wp) :: rn_zqt ! z(q,t) : height of humidity and temperature measurements 114 REAL(wp) :: rn_zu ! z(u) : height of wind measurements 115 !!gm ref namelist initialize it so remove the setting to false below 116 LOGICAL :: ln_Cd_L12 = .FALSE. ! Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2012) 117 LOGICAL :: ln_Cd_L15 = .FALSE. ! Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2015) 102 LOGICAL :: ln_Cd_L12 ! ice-atm drag = F( ice concentration ) (Lupkes et al. JGR2012) 103 LOGICAL :: ln_Cd_L15 ! ice-atm drag = F( ice concentration, atmospheric stability ) (Lupkes et al. JGR2015) 118 104 ! 119 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: Cd_atm ! transfer coefficient for momentum (tau) 120 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: Ch_atm ! transfer coefficient for sensible heat (Q_sens) 121 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: Ce_atm ! tansfert coefficient for evaporation (Q_lat) 122 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_zu ! air temperature at wind speed height (needed by Lupkes 2015 bulk scheme) 123 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: q_zu ! air spec. hum. at wind speed height (needed by Lupkes 2015 bulk scheme) 124 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: cdn_oce, chn_oce, cen_oce ! needed by Lupkes 2015 bulk scheme 105 LOGICAL :: ln_crt_fbk ! Add surface current feedback to the wind stress computation (Renault et al. 2020) 106 REAL(wp) :: rn_stau_a ! Alpha and Beta coefficients of Renault et al. 2020, eq. 10: Stau = Alpha * Wnd + Beta 107 REAL(wp) :: rn_stau_b ! 108 ! 109 REAL(wp) :: rn_pfac ! multiplication factor for precipitation 110 REAL(wp), PUBLIC :: rn_efac ! multiplication factor for evaporation 111 REAL(wp) :: rn_zqt ! z(q,t) : height of humidity and temperature measurements 112 REAL(wp) :: rn_zu ! z(u) : height of wind measurements 113 ! 114 REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: Cdn_oce, Chn_oce, Cen_oce ! neutral coeffs over ocean (L15 bulk scheme and ABL) 115 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: Cd_ice , Ch_ice , Ce_ice ! transfert coefficients over ice 116 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: t_zu, q_zu ! air temp. and spec. hum. at wind speed height (L15 bulk scheme) 117 118 LOGICAL :: ln_skin_cs ! use the cool-skin (only available in ECMWF and COARE algorithms) !LB 119 LOGICAL :: ln_skin_wl ! use the warm-layer parameterization (only available in ECMWF and COARE algorithms) !LB 120 LOGICAL :: ln_humi_sph ! humidity read in files ("sn_humi") is specific humidity [kg/kg] if .true. !LB 121 LOGICAL :: ln_humi_dpt ! humidity read in files ("sn_humi") is dew-point temperature [K] if .true. !LB 122 LOGICAL :: ln_humi_rlh ! humidity read in files ("sn_humi") is relative humidity [%] if .true. !LB 123 LOGICAL :: ln_tpot !!GS: flag to compute or not potential temperature 124 ! 125 INTEGER :: nhumi ! choice of the bulk algorithm 126 ! ! associated indices: 127 INTEGER, PARAMETER :: np_humi_sph = 1 128 INTEGER, PARAMETER :: np_humi_dpt = 2 129 INTEGER, PARAMETER :: np_humi_rlh = 3 125 130 126 131 INTEGER :: nblk ! choice of the bulk algorithm … … 128 133 INTEGER, PARAMETER :: np_NCAR = 1 ! "NCAR" algorithm (Large and Yeager 2008) 129 134 INTEGER, PARAMETER :: np_COARE_3p0 = 2 ! "COARE 3.0" algorithm (Fairall et al. 2003) 130 INTEGER, PARAMETER :: np_COARE_3p 5 = 3 ! "COARE 3.5" algorithm (Edson et al. 2013)131 INTEGER, PARAMETER :: np_ECMWF = 4 ! "ECMWF" algorithm (IFS cycle 31)135 INTEGER, PARAMETER :: np_COARE_3p6 = 3 ! "COARE 3.6" algorithm (Edson et al. 2013) 136 INTEGER, PARAMETER :: np_ECMWF = 4 ! "ECMWF" algorithm (IFS cycle 45r1) 132 137 133 138 !! * Substitutions 134 # include " vectopt_loop_substitute.h90"139 # include "do_loop_substitute.h90" 135 140 !!---------------------------------------------------------------------- 136 141 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 144 149 !! *** ROUTINE sbc_blk_alloc *** 145 150 !!------------------------------------------------------------------- 146 ALLOCATE( Cd_atm (jpi,jpj), Ch_atm (jpi,jpj), Ce_atm (jpi,jpj), t_zu(jpi,jpj), q_zu(jpi,jpj), & 147 & cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), STAT=sbc_blk_alloc ) 151 ALLOCATE( t_zu(jpi,jpj) , q_zu(jpi,jpj) , & 152 & Cdn_oce(jpi,jpj), Chn_oce(jpi,jpj), Cen_oce(jpi,jpj), & 153 & Cd_ice (jpi,jpj), Ch_ice (jpi,jpj), Ce_ice (jpi,jpj), STAT=sbc_blk_alloc ) 148 154 ! 149 155 CALL mpp_sum ( 'sbcblk', sbc_blk_alloc ) … … 158 164 !! ** Purpose : choose and initialize a bulk formulae formulation 159 165 !! 160 !! ** Method : 166 !! ** Method : 161 167 !! 162 168 !!---------------------------------------------------------------------- 163 INTEGER :: ifpr, jfld! dummy loop indice and argument169 INTEGER :: jfpr ! dummy loop indice and argument 164 170 INTEGER :: ios, ierror, ioptio ! Local integer 165 171 !! 166 172 CHARACTER(len=100) :: cn_dir ! Root directory for location of atmospheric forcing files 167 173 TYPE(FLD_N), DIMENSION(jpfld) :: slf_i ! array of namelist informations on the fields to read 168 TYPE(FLD_N) :: sn_wndi, sn_wndj, sn_humi, sn_qsr ! informations about the fields to be read 169 TYPE(FLD_N) :: sn_qlw , sn_tair, sn_prec, sn_snow ! " " 170 TYPE(FLD_N) :: sn_slp , sn_tdif ! " " 174 TYPE(FLD_N) :: sn_wndi, sn_wndj , sn_humi, sn_qsr ! informations about the fields to be read 175 TYPE(FLD_N) :: sn_qlw , sn_tair , sn_prec, sn_snow ! " " 176 TYPE(FLD_N) :: sn_slp , sn_uoatm, sn_voatm ! " " 177 TYPE(FLD_N) :: sn_hpgi, sn_hpgj ! " " 178 INTEGER :: ipka ! number of levels in the atmospheric variable 171 179 NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw , & ! input fields 172 & sn_tair, sn_prec, sn_snow, sn_slp, sn_tdif, & 173 & ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF, & ! bulk algorithm 174 & cn_dir , ln_taudif, rn_zqt, rn_zu, & 175 & rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15 180 & sn_tair, sn_prec, sn_snow, sn_slp, sn_uoatm, sn_voatm, & 181 & sn_hpgi, sn_hpgj, & 182 & ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF, & ! bulk algorithm 183 & cn_dir , rn_zqt, rn_zu, & 184 & rn_pfac, rn_efac, ln_Cd_L12, ln_Cd_L15, ln_tpot, & 185 & ln_crt_fbk, rn_stau_a, rn_stau_b, & ! current feedback 186 & ln_skin_cs, ln_skin_wl, ln_humi_sph, ln_humi_dpt, ln_humi_rlh ! cool-skin / warm-layer !LB 176 187 !!--------------------------------------------------------------------- 177 188 ! … … 179 190 IF( sbc_blk_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' ) 180 191 ! 181 ! !** read bulk namelist 182 REWIND( numnam_ref ) !* Namelist namsbc_blk in reference namelist : bulk parameters 192 ! !** read bulk namelist 183 193 READ ( numnam_ref, namsbc_blk, IOSTAT = ios, ERR = 901) 184 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_blk in reference namelist', lwp ) 185 ! 186 REWIND( numnam_cfg ) !* Namelist namsbc_blk in configuration namelist : bulk parameters 194 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_blk in reference namelist' ) 195 ! 187 196 READ ( numnam_cfg, namsbc_blk, IOSTAT = ios, ERR = 902 ) 188 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namsbc_blk in configuration namelist' , lwp)197 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namsbc_blk in configuration namelist' ) 189 198 ! 190 199 IF(lwm) WRITE( numond, namsbc_blk ) … … 192 201 ! !** initialization of the chosen bulk formulae (+ check) 193 202 ! !* select the bulk chosen in the namelist and check the choice 194 ioptio = 0 195 IF( ln_NCAR ) THEN ; nblk = np_NCAR ; ioptio = ioptio + 1 ; ENDIF 196 IF( ln_COARE_3p0 ) THEN ; nblk = np_COARE_3p0 ; ioptio = ioptio + 1 ; ENDIF 197 IF( ln_COARE_3p5 ) THEN ; nblk = np_COARE_3p5 ; ioptio = ioptio + 1 ; ENDIF 198 IF( ln_ECMWF ) THEN ; nblk = np_ECMWF ; ioptio = ioptio + 1 ; ENDIF 199 ! 203 ioptio = 0 204 IF( ln_NCAR ) THEN 205 nblk = np_NCAR ; ioptio = ioptio + 1 206 ENDIF 207 IF( ln_COARE_3p0 ) THEN 208 nblk = np_COARE_3p0 ; ioptio = ioptio + 1 209 ENDIF 210 IF( ln_COARE_3p6 ) THEN 211 nblk = np_COARE_3p6 ; ioptio = ioptio + 1 212 ENDIF 213 IF( ln_ECMWF ) THEN 214 nblk = np_ECMWF ; ioptio = ioptio + 1 215 ENDIF 200 216 IF( ioptio /= 1 ) CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' ) 217 218 ! !** initialization of the cool-skin / warm-layer parametrization 219 IF( ln_skin_cs .OR. ln_skin_wl ) THEN 220 !! Some namelist sanity tests: 221 IF( ln_NCAR ) & 222 & CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with NCAR algorithm' ) 223 IF( nn_fsbc /= 1 ) & 224 & CALL ctl_stop( 'sbc_blk_init: Please set "nn_fsbc" to 1 when using cool-skin/warm-layer param.') 225 END IF 226 227 IF( ln_skin_wl ) THEN 228 !! Check if the frequency of downwelling solar flux input makes sense and if ln_dm2dc=T if it is daily! 229 IF( (sn_qsr%freqh < 0.).OR.(sn_qsr%freqh > 24.) ) & 230 & CALL ctl_stop( 'sbc_blk_init: Warm-layer param. (ln_skin_wl) not compatible with freq. of solar flux > daily' ) 231 IF( (sn_qsr%freqh == 24.).AND.(.NOT. ln_dm2dc) ) & 232 & CALL ctl_stop( 'sbc_blk_init: Please set ln_dm2dc=T for warm-layer param. (ln_skin_wl) to work properly' ) 233 END IF 234 235 ioptio = 0 236 IF( ln_humi_sph ) THEN 237 nhumi = np_humi_sph ; ioptio = ioptio + 1 238 ENDIF 239 IF( ln_humi_dpt ) THEN 240 nhumi = np_humi_dpt ; ioptio = ioptio + 1 241 ENDIF 242 IF( ln_humi_rlh ) THEN 243 nhumi = np_humi_rlh ; ioptio = ioptio + 1 244 ENDIF 245 IF( ioptio /= 1 ) CALL ctl_stop( 'sbc_blk_init: Choose one and only one type of air humidity' ) 201 246 ! 202 247 IF( ln_dm2dc ) THEN !* check: diurnal cycle on Qsr 203 IF( sn_qsr% nfreqh /= 24) CALL ctl_stop( 'sbc_blk_init: ln_dm2dc=T only with daily short-wave input' )204 IF( sn_qsr%ln_tint ) THEN 248 IF( sn_qsr%freqh /= 24. ) CALL ctl_stop( 'sbc_blk_init: ln_dm2dc=T only with daily short-wave input' ) 249 IF( sn_qsr%ln_tint ) THEN 205 250 CALL ctl_warn( 'sbc_blk_init: ln_dm2dc=T daily qsr time interpolation done by sbcdcy module', & 206 251 & ' ==> We force time interpolation = .false. for qsr' ) … … 210 255 ! !* set the bulk structure 211 256 ! !- store namelist information in an array 212 slf_i(jp_wndi) = sn_wndi ; slf_i(jp_wndj) = sn_wndj 213 slf_i(jp_qsr ) = sn_qsr ; slf_i(jp_qlw ) = sn_qlw 214 slf_i(jp_tair) = sn_tair ; slf_i(jp_humi) = sn_humi 215 slf_i(jp_prec) = sn_prec ; slf_i(jp_snow) = sn_snow 216 slf_i(jp_slp) = sn_slp ; slf_i(jp_tdif) = sn_tdif 217 ! 218 lhftau = ln_taudif !- add an extra field if HF stress is used 219 jfld = jpfld - COUNT( (/.NOT.lhftau/) ) 257 ! 258 slf_i(jp_wndi ) = sn_wndi ; slf_i(jp_wndj ) = sn_wndj 259 slf_i(jp_qsr ) = sn_qsr ; slf_i(jp_qlw ) = sn_qlw 260 slf_i(jp_tair ) = sn_tair ; slf_i(jp_humi ) = sn_humi 261 slf_i(jp_prec ) = sn_prec ; slf_i(jp_snow ) = sn_snow 262 slf_i(jp_slp ) = sn_slp 263 slf_i(jp_uoatm) = sn_uoatm ; slf_i(jp_voatm) = sn_voatm 264 slf_i(jp_hpgi ) = sn_hpgi ; slf_i(jp_hpgj ) = sn_hpgj 265 ! 266 IF( .NOT. ln_abl ) THEN ! force to not use jp_hpgi and jp_hpgj, should already be done in namelist_* but we never know... 267 slf_i(jp_hpgi)%clname = 'NOT USED' 268 slf_i(jp_hpgj)%clname = 'NOT USED' 269 ENDIF 220 270 ! 221 271 ! !- allocate the bulk structure 222 ALLOCATE( sf(j fld), STAT=ierror )272 ALLOCATE( sf(jpfld), STAT=ierror ) 223 273 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_blk_init: unable to allocate sf structure' ) 224 DO ifpr= 1, jfld 225 ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) ) 226 IF( slf_i(ifpr)%ln_tint ) ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) ) 227 IF( slf_i(ifpr)%nfreqh > 0. .AND. MOD( 3600. * slf_i(ifpr)%nfreqh , REAL(nn_fsbc) * rdt) /= 0. ) & 228 & CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rdt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.', & 229 & ' This is not ideal. You should consider changing either rdt or nn_fsbc value...' ) 230 231 END DO 274 ! 232 275 ! !- fill the bulk structure with namelist informations 233 276 CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' ) 234 277 ! 235 IF ( ln_wave ) THEN 236 !Activated wave module but neither drag nor stokes drift activated 237 IF ( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) ) THEN 278 DO jfpr= 1, jpfld 279 ! 280 IF( ln_abl .AND. & 281 & ( jfpr == jp_wndi .OR. jfpr == jp_wndj .OR. jfpr == jp_humi .OR. & 282 & jfpr == jp_hpgi .OR. jfpr == jp_hpgj .OR. jfpr == jp_tair ) ) THEN 283 ipka = jpka ! ABL: some fields are 3D input 284 ELSE 285 ipka = 1 286 ENDIF 287 ! 288 ALLOCATE( sf(jfpr)%fnow(jpi,jpj,ipka) ) 289 ! 290 IF( TRIM(sf(jfpr)%clrootname) == 'NOT USED' ) THEN !-- not used field --! (only now allocated and set to default) 291 IF( jfpr == jp_slp ) THEN 292 sf(jfpr)%fnow(:,:,1:ipka) = 101325._wp ! use standard pressure in Pa 293 ELSEIF( jfpr == jp_prec .OR. jfpr == jp_snow .OR. jfpr == jp_uoatm .OR. jfpr == jp_voatm ) THEN 294 sf(jfpr)%fnow(:,:,1:ipka) = 0._wp ! no precip or no snow or no surface currents 295 ELSEIF( ( jfpr == jp_hpgi .OR. jfpr == jp_hpgj ) .AND. .NOT. ln_abl ) THEN 296 DEALLOCATE( sf(jfpr)%fnow ) ! deallocate as not used in this case 297 ELSE 298 WRITE(ctmp1,*) 'sbc_blk_init: no default value defined for field number', jfpr 299 CALL ctl_stop( ctmp1 ) 300 ENDIF 301 ELSE !-- used field --! 302 IF( sf(jfpr)%ln_tint ) ALLOCATE( sf(jfpr)%fdta(jpi,jpj,ipka,2) ) ! allocate array for temporal interpolation 303 ! 304 IF( sf(jfpr)%freqh > 0. .AND. MOD( NINT(3600. * sf(jfpr)%freqh), nn_fsbc * NINT(rn_Dt) ) /= 0 ) & 305 & CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rn_Dt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.', & 306 & ' This is not ideal. You should consider changing either rn_Dt or nn_fsbc value...' ) 307 ENDIF 308 END DO 309 ! 310 IF( ln_wave ) THEN 311 !Activated wave module but neither drag nor stokes drift activated 312 IF( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) ) THEN 238 313 CALL ctl_stop( 'STOP', 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauwoc=F, ln_stcor=F' ) 239 !drag coefficient read from wave model definable only with mfs bulk formulae and core240 ELSEIF (ln_cdgw .AND. .NOT. ln_NCAR ) THEN241 242 ELSEIF 243 314 !drag coefficient read from wave model definable only with mfs bulk formulae and core 315 ELSEIF(ln_cdgw .AND. .NOT. ln_NCAR ) THEN 316 CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR and CORE bulk formulae') 317 ELSEIF(ln_stcor .AND. .NOT. ln_sdw) THEN 318 CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T') 244 319 ENDIF 245 320 ELSE 246 IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) & 247 & CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ', & 248 & 'with drag coefficient (ln_cdgw =T) ' , & 249 & 'or Stokes Drift (ln_sdw=T) ' , & 250 & 'or ocean stress modification due to waves (ln_tauwoc=T) ', & 251 & 'or Stokes-Coriolis term (ln_stcori=T)' ) 252 ENDIF 253 ! 254 ! 321 IF( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) & 322 & CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ', & 323 & 'with drag coefficient (ln_cdgw =T) ' , & 324 & 'or Stokes Drift (ln_sdw=T) ' , & 325 & 'or ocean stress modification due to waves (ln_tauwoc=T) ', & 326 & 'or Stokes-Coriolis term (ln_stcori=T)' ) 327 ENDIF 328 ! 329 IF( ln_abl ) THEN ! ABL: read 3D fields for wind, temperature, humidity and pressure gradient 330 rn_zqt = ght_abl(2) ! set the bulk altitude to ABL first level 331 rn_zu = ght_abl(2) 332 IF(lwp) WRITE(numout,*) 333 IF(lwp) WRITE(numout,*) ' ABL formulation: overwrite rn_zqt & rn_zu with ABL first level altitude' 334 ENDIF 335 ! 336 ! set transfer coefficients to default sea-ice values 337 Cd_ice(:,:) = rCd_ice 338 Ch_ice(:,:) = rCd_ice 339 Ce_ice(:,:) = rCd_ice 340 ! 255 341 IF(lwp) THEN !** Control print 256 342 ! 257 WRITE(numout,*) !* namelist 343 WRITE(numout,*) !* namelist 258 344 WRITE(numout,*) ' Namelist namsbc_blk (other than data information):' 259 345 WRITE(numout,*) ' "NCAR" algorithm (Large and Yeager 2008) ln_NCAR = ', ln_NCAR 260 346 WRITE(numout,*) ' "COARE 3.0" algorithm (Fairall et al. 2003) ln_COARE_3p0 = ', ln_COARE_3p0 261 WRITE(numout,*) ' "COARE 3.5" algorithm (Edson et al. 2013) ln_COARE_3p5 = ', ln_COARE_3p0 262 WRITE(numout,*) ' "ECMWF" algorithm (IFS cycle 31) ln_ECMWF = ', ln_ECMWF 263 WRITE(numout,*) ' add High freq.contribution to the stress module ln_taudif = ', ln_taudif 347 WRITE(numout,*) ' "COARE 3.6" algorithm (Fairall 2018 + Edson al 2013)ln_COARE_3p6 = ', ln_COARE_3p6 348 WRITE(numout,*) ' "ECMWF" algorithm (IFS cycle 45r1) ln_ECMWF = ', ln_ECMWF 264 349 WRITE(numout,*) ' Air temperature and humidity reference height (m) rn_zqt = ', rn_zqt 265 350 WRITE(numout,*) ' Wind vector reference height (m) rn_zu = ', rn_zu 266 351 WRITE(numout,*) ' factor applied on precipitation (total & snow) rn_pfac = ', rn_pfac 267 352 WRITE(numout,*) ' factor applied on evaporation rn_efac = ', rn_efac 268 WRITE(numout,*) ' factor applied on ocean/ice velocity rn_vfac = ', rn_vfac269 353 WRITE(numout,*) ' (form absolute (=0) to relative winds(=1))' 270 354 WRITE(numout,*) ' use ice-atm drag from Lupkes2012 ln_Cd_L12 = ', ln_Cd_L12 271 355 WRITE(numout,*) ' use ice-atm drag from Lupkes2015 ln_Cd_L15 = ', ln_Cd_L15 356 WRITE(numout,*) ' use surface current feedback on wind stress ln_crt_fbk = ', ln_crt_fbk 357 IF(ln_crt_fbk) THEN 358 WRITE(numout,*) ' Renault et al. 2020, eq. 10: Stau = Alpha * Wnd + Beta' 359 WRITE(numout,*) ' Alpha rn_stau_a = ', rn_stau_a 360 WRITE(numout,*) ' Beta rn_stau_b = ', rn_stau_b 361 ENDIF 272 362 ! 273 363 WRITE(numout,*) … … 275 365 CASE( np_NCAR ) ; WRITE(numout,*) ' ==>>> "NCAR" algorithm (Large and Yeager 2008)' 276 366 CASE( np_COARE_3p0 ) ; WRITE(numout,*) ' ==>>> "COARE 3.0" algorithm (Fairall et al. 2003)' 277 CASE( np_COARE_3p5 ) ; WRITE(numout,*) ' ==>>> "COARE 3.5" algorithm (Edson et al. 2013)' 278 CASE( np_ECMWF ) ; WRITE(numout,*) ' ==>>> "ECMWF" algorithm (IFS cycle 31)' 367 CASE( np_COARE_3p6 ) ; WRITE(numout,*) ' ==>>> "COARE 3.6" algorithm (Fairall 2018+Edson et al. 2013)' 368 CASE( np_ECMWF ) ; WRITE(numout,*) ' ==>>> "ECMWF" algorithm (IFS cycle 45r1)' 369 END SELECT 370 ! 371 WRITE(numout,*) 372 WRITE(numout,*) ' use cool-skin parameterization (SSST) ln_skin_cs = ', ln_skin_cs 373 WRITE(numout,*) ' use warm-layer parameterization (SSST) ln_skin_wl = ', ln_skin_wl 374 ! 375 WRITE(numout,*) 376 SELECT CASE( nhumi ) !* Print the choice of air humidity 377 CASE( np_humi_sph ) ; WRITE(numout,*) ' ==>>> air humidity is SPECIFIC HUMIDITY [kg/kg]' 378 CASE( np_humi_dpt ) ; WRITE(numout,*) ' ==>>> air humidity is DEW-POINT TEMPERATURE [K]' 379 CASE( np_humi_rlh ) ; WRITE(numout,*) ' ==>>> air humidity is RELATIVE HUMIDITY [%]' 279 380 END SELECT 280 381 ! … … 297 398 !! (momentum, heat, freshwater and runoff) 298 399 !! 299 !! ** Method : (1) READ each fluxes in NetCDF files: 300 !! the 10m wind velocity (i-component) (m/s) at T-point 301 !! the 10m wind velocity (j-component) (m/s) at T-point 302 !! the 10m or 2m specific humidity ( % ) 303 !! the solar heat (W/m2) 304 !! the Long wave (W/m2) 305 !! the 10m or 2m air temperature (Kelvin) 306 !! the total precipitation (rain+snow) (Kg/m2/s) 307 !! the snow (solid prcipitation) (kg/m2/s) 308 !! the tau diff associated to HF tau (N/m2) at T-point (ln_taudif=T) 309 !! (2) CALL blk_oce 400 !! ** Method : 401 !! (1) READ each fluxes in NetCDF files: 402 !! the wind velocity (i-component) at z=rn_zu (m/s) at T-point 403 !! the wind velocity (j-component) at z=rn_zu (m/s) at T-point 404 !! the specific humidity at z=rn_zqt (kg/kg) 405 !! the air temperature at z=rn_zqt (Kelvin) 406 !! the solar heat (W/m2) 407 !! the Long wave (W/m2) 408 !! the total precipitation (rain+snow) (Kg/m2/s) 409 !! the snow (solid precipitation) (kg/m2/s) 410 !! ABL dynamical forcing (i/j-components of either hpg or geostrophic winds) 411 !! (2) CALL blk_oce_1 and blk_oce_2 310 412 !! 311 413 !! C A U T I O N : never mask the surface stress fields … … 324 426 !!---------------------------------------------------------------------- 325 427 INTEGER, INTENT(in) :: kt ! ocean time step 326 !!--------------------------------------------------------------------- 428 !!---------------------------------------------------------------------- 429 REAL(wp), DIMENSION(jpi,jpj) :: zssq, zcd_du, zsen, zevp 430 REAL(wp) :: ztmp 431 !!---------------------------------------------------------------------- 327 432 ! 328 433 CALL fld_read( kt, nn_fsbc, sf ) ! input fields provided at the current time-step 329 ! 434 435 ! Sanity/consistence test on humidity at first time step to detect potential screw-up: 436 IF( kt == nit000 ) THEN 437 IF(lwp) WRITE(numout,*) '' 438 #if defined key_agrif 439 IF(lwp) WRITE(numout,*) ' === AGRIF => Sanity/consistence test on air humidity SKIPPED! :( ===' 440 #else 441 ztmp = SUM(tmask(:,:,1)) ! number of ocean points on local proc domain 442 IF( ztmp > 8._wp ) THEN ! test only on proc domains with at least 8 ocean points! 443 ztmp = SUM(sf(jp_humi)%fnow(:,:,1)*tmask(:,:,1))/ztmp ! mean humidity over ocean on proc 444 SELECT CASE( nhumi ) 445 CASE( np_humi_sph ) ! specific humidity => expect: 0. <= something < 0.065 [kg/kg] (0.061 is saturation at 45degC !!!) 446 IF( (ztmp < 0._wp) .OR. (ztmp > 0.065) ) ztmp = -1._wp 447 CASE( np_humi_dpt ) ! dew-point temperature => expect: 110. <= something < 320. [K] 448 IF( (ztmp < 110._wp).OR.(ztmp > 320._wp) ) ztmp = -1._wp 449 CASE( np_humi_rlh ) ! relative humidity => expect: 0. <= something < 100. [%] 450 IF( (ztmp < 0._wp) .OR.(ztmp > 100._wp) ) ztmp = -1._wp 451 END SELECT 452 IF(ztmp < 0._wp) THEN 453 IF (lwp) WRITE(numout,'(" Mean humidity value found on proc #",i6.6," is: ",f10.5)') narea, ztmp 454 CALL ctl_stop( 'STOP', 'Something is wrong with air humidity!!!', & 455 & ' ==> check the unit in your input files' , & 456 & ' ==> check consistence of namelist choice: specific? relative? dew-point?', & 457 & ' ==> ln_humi_sph -> [kg/kg] | ln_humi_rlh -> [%] | ln_humi_dpt -> [K] !!!' ) 458 END IF 459 END IF 460 IF(lwp) WRITE(numout,*) ' === Sanity/consistence test on air humidity sucessfuly passed! ===' 461 #endif 462 IF(lwp) WRITE(numout,*) '' 463 END IF !IF( kt == nit000 ) 330 464 ! ! compute the surface ocean fluxes using bulk formulea 331 IF( MOD( kt - 1, nn_fsbc ) == 0 ) CALL blk_oce( kt, sf, sst_m, ssu_m, ssv_m ) 332 465 IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 466 CALL blk_oce_1( kt, sf(jp_wndi )%fnow(:,:,1), sf(jp_wndj )%fnow(:,:,1), & ! <<= in 467 & sf(jp_tair )%fnow(:,:,1), sf(jp_humi )%fnow(:,:,1), & ! <<= in 468 & sf(jp_slp )%fnow(:,:,1), sst_m, ssu_m, ssv_m, & ! <<= in 469 & sf(jp_uoatm)%fnow(:,:,1), sf(jp_voatm)%fnow(:,:,1), & ! <<= in 470 & sf(jp_qsr )%fnow(:,:,1), sf(jp_qlw )%fnow(:,:,1), & ! <<= in (wl/cs) 471 & tsk_m, zssq, zcd_du, zsen, zevp ) ! =>> out 472 473 CALL blk_oce_2( sf(jp_tair )%fnow(:,:,1), sf(jp_qsr )%fnow(:,:,1), & ! <<= in 474 & sf(jp_qlw )%fnow(:,:,1), sf(jp_prec )%fnow(:,:,1), & ! <<= in 475 & sf(jp_snow )%fnow(:,:,1), tsk_m, & ! <<= in 476 & zsen, zevp ) ! <=> in out 477 ENDIF 478 ! 333 479 #if defined key_cice 334 480 IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 335 481 qlw_ice(:,:,1) = sf(jp_qlw )%fnow(:,:,1) 336 IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) 337 ELSE ; qsr_ice(:,:,1) = sf(jp_qsr)%fnow(:,:,1) 338 ENDIF 482 IF( ln_dm2dc ) THEN 483 qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) 484 ELSE 485 qsr_ice(:,:,1) = sf(jp_qsr)%fnow(:,:,1) 486 ENDIF 339 487 tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) 340 qatm_ice(:,:) = sf(jp_humi)%fnow(:,:,1) 488 489 SELECT CASE( nhumi ) 490 CASE( np_humi_sph ) 491 qatm_ice(:,:) = sf(jp_humi)%fnow(:,:,1) 492 CASE( np_humi_dpt ) 493 qatm_ice(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 494 CASE( np_humi_rlh ) 495 qatm_ice(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) !LB: 0.01 => RH is % percent in file 496 END SELECT 497 341 498 tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac 342 499 sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac … … 349 506 350 507 351 SUBROUTINE blk_oce( kt, sf, pst, pu, pv ) 352 !!--------------------------------------------------------------------- 353 !! *** ROUTINE blk_oce *** 354 !! 355 !! ** Purpose : provide the momentum, heat and freshwater fluxes at 356 !! the ocean surface at each time step 357 !! 358 !! ** Method : bulk formulea for the ocean using atmospheric 359 !! fields read in sbc_read 508 SUBROUTINE blk_oce_1( kt, pwndi, pwndj, ptair, phumi, & ! inp 509 & pslp , pst , pu , pv, & ! inp 510 & puatm, pvatm, pqsr , pqlw , & ! inp 511 & ptsk , pssq , pcd_du, psen, pevp ) ! out 512 !!--------------------------------------------------------------------- 513 !! *** ROUTINE blk_oce_1 *** 514 !! 515 !! ** Purpose : if ln_blk=T, computes surface momentum, heat and freshwater fluxes 516 !! if ln_abl=T, computes Cd x |U|, Ch x |U|, Ce x |U| for ABL integration 517 !! 518 !! ** Method : bulk formulae using atmospheric fields from : 519 !! if ln_blk=T, atmospheric fields read in sbc_read 520 !! if ln_abl=T, the ABL model at previous time-step 521 !! 522 !! ** Outputs : - pssq : surface humidity used to compute latent heat flux (kg/kg) 523 !! - pcd_du : Cd x |dU| at T-points (m/s) 524 !! - psen : Ch x |dU| at T-points (m/s) 525 !! - pevp : Ce x |dU| at T-points (m/s) 526 !!--------------------------------------------------------------------- 527 INTEGER , INTENT(in ) :: kt ! time step index 528 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pwndi ! atmospheric wind at U-point [m/s] 529 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pwndj ! atmospheric wind at V-point [m/s] 530 REAL(wp), INTENT(in ), DIMENSION(:,:) :: phumi ! specific humidity at T-points [kg/kg] 531 REAL(wp), INTENT(in ), DIMENSION(:,:) :: ptair ! potential temperature at T-points [Kelvin] 532 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pslp ! sea-level pressure [Pa] 533 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pst ! surface temperature [Celsius] 534 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pu ! surface current at U-point (i-component) [m/s] 535 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pv ! surface current at V-point (j-component) [m/s] 536 REAL(wp), INTENT(in ), DIMENSION(:,:) :: puatm ! surface current seen by the atm at T-point (i-component) [m/s] 537 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pvatm ! surface current seen by the atm at T-point (j-component) [m/s] 538 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pqsr ! 539 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pqlw ! 540 REAL(wp), INTENT( out), DIMENSION(:,:) :: ptsk ! skin temp. (or SST if CS & WL not used) [Celsius] 541 REAL(wp), INTENT( out), DIMENSION(:,:) :: pssq ! specific humidity at pst [kg/kg] 542 REAL(wp), INTENT( out), DIMENSION(:,:) :: pcd_du ! Cd x |dU| at T-points [m/s] 543 REAL(wp), INTENT( out), DIMENSION(:,:) :: psen ! Ch x |dU| at T-points [m/s] 544 REAL(wp), INTENT( out), DIMENSION(:,:) :: pevp ! Ce x |dU| at T-points [m/s] 545 ! 546 INTEGER :: ji, jj ! dummy loop indices 547 REAL(wp) :: zztmp ! local variable 548 REAL(wp) :: zstmax, zstau 549 #if defined key_cyclone 550 REAL(wp), DIMENSION(jpi,jpj) :: zwnd_i, zwnd_j ! wind speed components at T-point 551 #endif 552 REAL(wp), DIMENSION(jpi,jpj) :: ztau_i, ztau_j ! wind stress components at T-point 553 REAL(wp), DIMENSION(jpi,jpj) :: zU_zu ! bulk wind speed at height zu [m/s] 554 REAL(wp), DIMENSION(jpi,jpj) :: ztpot ! potential temperature of air at z=rn_zqt [K] 555 REAL(wp), DIMENSION(jpi,jpj) :: zqair ! specific humidity of air at z=rn_zqt [kg/kg] 556 REAL(wp), DIMENSION(jpi,jpj) :: zcd_oce ! momentum transfert coefficient over ocean 557 REAL(wp), DIMENSION(jpi,jpj) :: zch_oce ! sensible heat transfert coefficient over ocean 558 REAL(wp), DIMENSION(jpi,jpj) :: zce_oce ! latent heat transfert coefficient over ocean 559 REAL(wp), DIMENSION(jpi,jpj) :: zqla ! latent heat flux 560 REAL(wp), DIMENSION(jpi,jpj) :: zztmp1, zztmp2 561 !!--------------------------------------------------------------------- 562 ! 563 ! local scalars ( place there for vector optimisation purposes) 564 ! ! Temporary conversion from Celcius to Kelvin (and set minimum value far above 0 K) 565 ptsk(:,:) = pst(:,:) + rt0 ! by default: skin temperature = "bulk SST" (will remain this way if NCAR algorithm used!) 566 567 ! ----------------------------------------------------------------------------- ! 568 ! 0 Wind components and module at T-point relative to the moving ocean ! 569 ! ----------------------------------------------------------------------------- ! 570 571 ! ... components ( U10m - U_oce ) at T-point (unmasked) 572 #if defined key_cyclone 573 zwnd_i(:,:) = 0._wp 574 zwnd_j(:,:) = 0._wp 575 CALL wnd_cyc( kt, zwnd_i, zwnd_j ) ! add analytical tropical cyclone (Vincent et al. JGR 2012) 576 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 577 zwnd_i(ji,jj) = pwndi(ji,jj) + zwnd_i(ji,jj) 578 zwnd_j(ji,jj) = pwndj(ji,jj) + zwnd_j(ji,jj) 579 ! ... scalar wind at T-point (not masked) 580 wndm(ji,jj) = SQRT( zwnd_i(ji,jj) * zwnd_i(ji,jj) + zwnd_j(ji,jj) * zwnd_j(ji,jj) ) 581 END_2D 582 #else 583 ! ... scalar wind module at T-point (not masked) 584 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 585 wndm(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) ) 586 END_2D 587 #endif 588 ! ----------------------------------------------------------------------------- ! 589 ! I Solar FLUX ! 590 ! ----------------------------------------------------------------------------- ! 591 592 ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle ! Short Wave 593 zztmp = 1. - albo 594 IF( ln_dm2dc ) THEN 595 qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1) 596 ELSE 597 qsr(:,:) = zztmp * sf(jp_qsr)%fnow(:,:,1) * tmask(:,:,1) 598 ENDIF 599 600 601 ! ----------------------------------------------------------------------------- ! 602 ! II Turbulent FLUXES ! 603 ! ----------------------------------------------------------------------------- ! 604 605 ! specific humidity at SST 606 pssq(:,:) = rdct_qsat_salt * q_sat( ptsk(:,:), pslp(:,:) ) 607 608 IF( ln_skin_cs .OR. ln_skin_wl ) THEN 609 !! Backup "bulk SST" and associated spec. hum. 610 zztmp1(:,:) = ptsk(:,:) 611 zztmp2(:,:) = pssq(:,:) 612 ENDIF 613 614 ! specific humidity of air at "rn_zqt" m above the sea 615 SELECT CASE( nhumi ) 616 CASE( np_humi_sph ) 617 zqair(:,:) = phumi(:,:) ! what we read in file is already a spec. humidity! 618 CASE( np_humi_dpt ) 619 !IF(lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of d_air and slp !' !LBrm 620 zqair(:,:) = q_sat( phumi(:,:), pslp(:,:) ) 621 CASE( np_humi_rlh ) 622 !IF(lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of RH, t_air and slp !' !LBrm 623 zqair(:,:) = q_air_rh( 0.01_wp*phumi(:,:), ptair(:,:), pslp(:,:) ) !LB: 0.01 => RH is % percent in file 624 END SELECT 625 ! 626 ! potential temperature of air at "rn_zqt" m above the sea 627 IF( ln_abl ) THEN 628 ztpot = ptair(:,:) 629 ELSE 630 ! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate 631 ! (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2 632 ! (since reanalysis products provide T at z, not theta !) 633 !#LB: because AGRIF hates functions that return something else than a scalar, need to 634 ! use scalar version of gamma_moist() ... 635 IF( ln_tpot ) THEN 636 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 637 ztpot(ji,jj) = ptair(ji,jj) + gamma_moist( ptair(ji,jj), zqair(ji,jj) ) * rn_zqt 638 END_2D 639 ELSE 640 ztpot = ptair(:,:) 641 ENDIF 642 ENDIF 643 644 !! Time to call the user-selected bulk parameterization for 645 !! == transfer coefficients ==! Cd, Ch, Ce at T-point, and more... 646 SELECT CASE( nblk ) 647 648 CASE( np_NCAR ) 649 CALL turb_ncar ( rn_zqt, rn_zu, ptsk, ztpot, pssq, zqair, wndm, & 650 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 651 652 CASE( np_COARE_3p0 ) 653 CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, ptsk, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 654 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 655 & Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 656 657 CASE( np_COARE_3p6 ) 658 CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, ptsk, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 659 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 660 & Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 661 662 CASE( np_ECMWF ) 663 CALL turb_ecmwf ( kt, rn_zqt, rn_zu, ptsk, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 664 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 665 & Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 666 667 CASE DEFAULT 668 CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 669 670 END SELECT 671 672 IF( iom_use('Cd_oce') ) CALL iom_put("Cd_oce", zcd_oce * tmask(:,:,1)) 673 IF( iom_use('Ce_oce') ) CALL iom_put("Ce_oce", zce_oce * tmask(:,:,1)) 674 IF( iom_use('Ch_oce') ) CALL iom_put("Ch_oce", zch_oce * tmask(:,:,1)) 675 !! LB: mainly here for debugging purpose: 676 IF( iom_use('theta_zt') ) CALL iom_put("theta_zt", (ztpot-rt0) * tmask(:,:,1)) ! potential temperature at z=zt 677 IF( iom_use('q_zt') ) CALL iom_put("q_zt", zqair * tmask(:,:,1)) ! specific humidity " 678 IF( iom_use('theta_zu') ) CALL iom_put("theta_zu", (t_zu -rt0) * tmask(:,:,1)) ! potential temperature at z=zu 679 IF( iom_use('q_zu') ) CALL iom_put("q_zu", q_zu * tmask(:,:,1)) ! specific humidity " 680 IF( iom_use('ssq') ) CALL iom_put("ssq", pssq * tmask(:,:,1)) ! saturation specific humidity at z=0 681 IF( iom_use('wspd_blk') ) CALL iom_put("wspd_blk", zU_zu * tmask(:,:,1)) ! bulk wind speed at z=zu 682 683 IF( ln_skin_cs .OR. ln_skin_wl ) THEN 684 !! ptsk and pssq have been updated!!! 685 !! 686 !! In the presence of sea-ice we forget about the cool-skin/warm-layer update of ptsk and pssq: 687 WHERE ( fr_i(:,:) > 0.001_wp ) 688 ! sea-ice present, we forget about the update, using what we backed up before call to turb_*() 689 ptsk(:,:) = zztmp1(:,:) 690 pssq(:,:) = zztmp2(:,:) 691 END WHERE 692 END IF 693 694 ! Turbulent fluxes over ocean => BULK_FORMULA @ sbcblk_phy.F90 695 ! ------------------------------------------------------------- 696 697 IF( ln_abl ) THEN !== ABL formulation ==! multiplication by rho_air and turbulent fluxes computation done in ablstp 698 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 699 zztmp = zU_zu(ji,jj) 700 wndm(ji,jj) = zztmp ! Store zU_zu in wndm to compute ustar2 in ablmod 701 pcd_du(ji,jj) = zztmp * zcd_oce(ji,jj) 702 psen(ji,jj) = zztmp * zch_oce(ji,jj) 703 pevp(ji,jj) = zztmp * zce_oce(ji,jj) 704 rhoa(ji,jj) = rho_air( ptair(ji,jj), phumi(ji,jj), pslp(ji,jj) ) 705 END_2D 706 ELSE !== BLK formulation ==! turbulent fluxes computation 707 CALL BULK_FORMULA( rn_zu, ptsk(:,:), pssq(:,:), t_zu(:,:), q_zu(:,:), & 708 & zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:), & 709 & wndm(:,:), zU_zu(:,:), pslp(:,:), & 710 & taum(:,:), psen(:,:), zqla(:,:), & 711 & pEvap=pevp(:,:), prhoa=rhoa(:,:), pfact_evap=rn_efac ) 712 713 zqla(:,:) = zqla(:,:) * tmask(:,:,1) 714 psen(:,:) = psen(:,:) * tmask(:,:,1) 715 taum(:,:) = taum(:,:) * tmask(:,:,1) 716 pevp(:,:) = pevp(:,:) * tmask(:,:,1) 717 718 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 719 IF( wndm(ji,jj) > 0._wp ) THEN 720 zztmp = taum(ji,jj) / wndm(ji,jj) 721 #if defined key_cyclone 722 ztau_i(ji,jj) = zztmp * zwnd_i(ji,jj) 723 ztau_j(ji,jj) = zztmp * zwnd_j(ji,jj) 724 #else 725 ztau_i(ji,jj) = zztmp * pwndi(ji,jj) 726 ztau_j(ji,jj) = zztmp * pwndj(ji,jj) 727 #endif 728 ELSE 729 ztau_i(ji,jj) = 0._wp 730 ztau_j(ji,jj) = 0._wp 731 ENDIF 732 END_2D 733 734 IF( ln_crt_fbk ) THEN ! aply eq. 10 and 11 of Renault et al. 2020 (doi: 10.1029/2019MS001715) 735 zstmax = MIN( rn_stau_a * 3._wp + rn_stau_b, 0._wp ) ! set the max value of Stau corresponding to a wind of 3 m/s (<0) 736 DO_2D( 0, 1, 0, 1 ) ! end at jpj and jpi, as ztau_j(ji,jj+1) ztau_i(ji+1,jj) used in the next loop 737 zstau = MIN( rn_stau_a * wndm(ji,jj) + rn_stau_b, zstmax ) ! stau (<0) must be smaller than zstmax 738 ztau_i(ji,jj) = ztau_i(ji,jj) + zstau * ( 0.5_wp * ( pu(ji-1,jj ) + pu(ji,jj) ) - puatm(ji,jj) ) 739 ztau_j(ji,jj) = ztau_j(ji,jj) + zstau * ( 0.5_wp * ( pv(ji ,jj-1) + pv(ji,jj) ) - pvatm(ji,jj) ) 740 taum(ji,jj) = SQRT( ztau_i(ji,jj) * ztau_i(ji,jj) + ztau_j(ji,jj) * ztau_j(ji,jj) ) 741 END_2D 742 ENDIF 743 744 ! ... utau, vtau at U- and V_points, resp. 745 ! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines 746 ! Note that coastal wind stress is not used in the code... so this extra care has no effect 747 DO_2D( 0, 0, 0, 0 ) ! start loop at 2, in case ln_crt_fbk = T 748 utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( ztau_i(ji,jj) + ztau_i(ji+1,jj ) ) & 749 & * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1)) 750 vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( ztau_j(ji,jj) + ztau_j(ji ,jj+1) ) & 751 & * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1)) 752 END_2D 753 754 IF( ln_crt_fbk ) THEN 755 CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1., taum, 'T', -1. ) 756 ELSE 757 CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. ) 758 ENDIF 759 760 CALL iom_put( "taum_oce", taum ) ! output wind stress module 761 762 IF(sn_cfctl%l_prtctl) THEN 763 CALL prt_ctl( tab2d_1=wndm , clinfo1=' blk_oce_1: wndm : ') 764 CALL prt_ctl( tab2d_1=utau , clinfo1=' blk_oce_1: utau : ', mask1=umask, & 765 & tab2d_2=vtau , clinfo2=' vtau : ', mask2=vmask ) 766 ENDIF 767 ! 768 ENDIF !IF( ln_abl ) 769 770 ptsk(:,:) = ( ptsk(:,:) - rt0 ) * tmask(:,:,1) ! Back to Celsius 771 772 IF( ln_skin_cs .OR. ln_skin_wl ) THEN 773 CALL iom_put( "t_skin" , ptsk ) ! T_skin in Celsius 774 CALL iom_put( "dt_skin" , ptsk - pst ) ! T_skin - SST temperature difference... 775 ENDIF 776 777 IF(sn_cfctl%l_prtctl) THEN 778 CALL prt_ctl( tab2d_1=pevp , clinfo1=' blk_oce_1: pevp : ' ) 779 CALL prt_ctl( tab2d_1=psen , clinfo1=' blk_oce_1: psen : ' ) 780 CALL prt_ctl( tab2d_1=pssq , clinfo1=' blk_oce_1: pssq : ' ) 781 ENDIF 782 ! 783 END SUBROUTINE blk_oce_1 784 785 786 SUBROUTINE blk_oce_2( ptair, pqsr, pqlw, pprec, & ! <<= in 787 & psnow, ptsk, psen, pevp ) ! <<= in 788 !!--------------------------------------------------------------------- 789 !! *** ROUTINE blk_oce_2 *** 790 !! 791 !! ** Purpose : finalize the momentum, heat and freshwater fluxes computation 792 !! at the ocean surface at each time step knowing Cd, Ch, Ce and 793 !! atmospheric variables (from ABL or external data) 360 794 !! 361 795 !! ** Outputs : - utau : i-component of the stress at U-point (N/m2) … … 366 800 !! - qns : Non Solar heat flux over the ocean (W/m2) 367 801 !! - emp : evaporation minus precipitation (kg/m2/s) 368 !! 369 !! ** Nota : sf has to be a dummy argument for AGRIF on NEC 370 !!--------------------------------------------------------------------- 371 INTEGER , INTENT(in ) :: kt ! time step index 372 TYPE(fld), INTENT(inout), DIMENSION(:) :: sf ! input data 373 REAL(wp) , INTENT(in) , DIMENSION(:,:) :: pst ! surface temperature [Celcius] 374 REAL(wp) , INTENT(in) , DIMENSION(:,:) :: pu ! surface current at U-point (i-component) [m/s] 375 REAL(wp) , INTENT(in) , DIMENSION(:,:) :: pv ! surface current at V-point (j-component) [m/s] 802 !!--------------------------------------------------------------------- 803 REAL(wp), INTENT(in), DIMENSION(:,:) :: ptair 804 REAL(wp), INTENT(in), DIMENSION(:,:) :: pqsr 805 REAL(wp), INTENT(in), DIMENSION(:,:) :: pqlw 806 REAL(wp), INTENT(in), DIMENSION(:,:) :: pprec 807 REAL(wp), INTENT(in), DIMENSION(:,:) :: psnow 808 REAL(wp), INTENT(in), DIMENSION(:,:) :: ptsk ! SKIN surface temperature [Celsius] 809 REAL(wp), INTENT(in), DIMENSION(:,:) :: psen 810 REAL(wp), INTENT(in), DIMENSION(:,:) :: pevp 376 811 ! 377 812 INTEGER :: ji, jj ! dummy loop indices 378 REAL(wp) :: zztmp ! local variable 379 REAL(wp), DIMENSION(jpi,jpj) :: zwnd_i, zwnd_j ! wind speed components at T-point 380 REAL(wp), DIMENSION(jpi,jpj) :: zsq ! specific humidity at pst 381 REAL(wp), DIMENSION(jpi,jpj) :: zqlw, zqsb ! long wave and sensible heat fluxes 382 REAL(wp), DIMENSION(jpi,jpj) :: zqla, zevap ! latent heat fluxes and evaporation 383 REAL(wp), DIMENSION(jpi,jpj) :: zst ! surface temperature in Kelvin 384 REAL(wp), DIMENSION(jpi,jpj) :: zU_zu ! bulk wind speed at height zu [m/s] 385 REAL(wp), DIMENSION(jpi,jpj) :: ztpot ! potential temperature of air at z=rn_zqt [K] 386 REAL(wp), DIMENSION(jpi,jpj) :: zrhoa ! density of air [kg/m^3] 813 REAL(wp) :: zztmp,zz1,zz2,zz3 ! local variable 814 REAL(wp), DIMENSION(jpi,jpj) :: ztskk ! skin temp. in Kelvin 815 REAL(wp), DIMENSION(jpi,jpj) :: zqlw ! long wave and sensible heat fluxes 816 REAL(wp), DIMENSION(jpi,jpj) :: zqla ! latent heat fluxes and evaporation 387 817 !!--------------------------------------------------------------------- 388 818 ! 389 819 ! local scalars ( place there for vector optimisation purposes) 390 zst(:,:) = pst(:,:) + rt0 ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K) 391 820 821 822 ztskk(:,:) = ptsk(:,:) + rt0 ! => ptsk in Kelvin rather than Celsius 823 392 824 ! ----------------------------------------------------------------------------- ! 393 ! 0 Wind components and module at T-point relative to the moving ocean!825 ! III Net longwave radiative FLUX ! 394 826 ! ----------------------------------------------------------------------------- ! 395 827 396 ! ... components ( U10m - U_oce ) at T-point (unmasked) 397 !!gm move zwnd_i (_j) set to zero inside the key_cyclone ??? 398 zwnd_i(:,:) = 0._wp 399 zwnd_j(:,:) = 0._wp 400 #if defined key_cyclone 401 CALL wnd_cyc( kt, zwnd_i, zwnd_j ) ! add analytical tropical cyclone (Vincent et al. JGR 2012) 402 DO jj = 2, jpjm1 403 DO ji = fs_2, fs_jpim1 ! vect. opt. 404 sf(jp_wndi)%fnow(ji,jj,1) = sf(jp_wndi)%fnow(ji,jj,1) + zwnd_i(ji,jj) 405 sf(jp_wndj)%fnow(ji,jj,1) = sf(jp_wndj)%fnow(ji,jj,1) + zwnd_j(ji,jj) 406 END DO 407 END DO 408 #endif 409 DO jj = 2, jpjm1 410 DO ji = fs_2, fs_jpim1 ! vect. opt. 411 zwnd_i(ji,jj) = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj ) + pu(ji,jj) ) ) 412 zwnd_j(ji,jj) = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji ,jj-1) + pv(ji,jj) ) ) 413 END DO 414 END DO 415 CALL lbc_lnk_multi( 'sbcblk', zwnd_i, 'T', -1., zwnd_j, 'T', -1. ) 416 ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked) 417 wndm(:,:) = SQRT( zwnd_i(:,:) * zwnd_i(:,:) & 418 & + zwnd_j(:,:) * zwnd_j(:,:) ) * tmask(:,:,1) 828 !! LB: now moved after Turbulent fluxes because must use the skin temperature rather that the SST 829 !! (ztskk is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.) 830 zqlw(:,:) = emiss_w * ( pqlw(:,:) - stefan*ztskk(:,:)*ztskk(:,:)*ztskk(:,:)*ztskk(:,:) ) * tmask(:,:,1) ! Net radiative longwave flux 831 832 ! Latent flux over ocean 833 ! ----------------------- 834 835 ! use scalar version of L_vap() for AGRIF compatibility 836 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 837 zqla(ji,jj) = - L_vap( ztskk(ji,jj) ) * pevp(ji,jj) ! Latent Heat flux !!GS: possibility to add a global qla to avoid recomputation after abl update 838 END_2D 839 840 IF(sn_cfctl%l_prtctl) THEN 841 CALL prt_ctl( tab2d_1=zqla , clinfo1=' blk_oce_2: zqla : ' ) 842 CALL prt_ctl( tab2d_1=zqlw , clinfo1=' blk_oce_2: zqlw : ', tab2d_2=qsr, clinfo2=' qsr : ' ) 843 844 ENDIF 419 845 420 846 ! ----------------------------------------------------------------------------- ! 421 ! I Radiative FLUXES!847 ! IV Total FLUXES ! 422 848 ! ----------------------------------------------------------------------------- ! 423 424 ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle ! Short Wave 425 zztmp = 1. - albo 426 IF( ln_dm2dc ) THEN ; qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1) 427 ELSE ; qsr(:,:) = zztmp * sf(jp_qsr)%fnow(:,:,1) * tmask(:,:,1) 428 ENDIF 429 430 zqlw(:,:) = ( sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:) ) * tmask(:,:,1) ! Long Wave 431 432 ! ----------------------------------------------------------------------------- ! 433 ! II Turbulent FLUXES ! 434 ! ----------------------------------------------------------------------------- ! 435 436 ! ... specific humidity at SST and IST tmask( 437 zsq(:,:) = 0.98 * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) ) 438 !! 439 !! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate 440 !! (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2 441 !! (since reanalysis products provide T at z, not theta !) 442 ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1) ) * rn_zqt 443 444 SELECT CASE( nblk ) !== transfer coefficients ==! Cd, Ch, Ce at T-point 445 ! 446 CASE( np_NCAR ) ; CALL turb_ncar ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! NCAR-COREv2 447 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 448 CASE( np_COARE_3p0 ) ; CALL turb_coare ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! COARE v3.0 449 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 450 CASE( np_COARE_3p5 ) ; CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! COARE v3.5 451 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 452 CASE( np_ECMWF ) ; CALL turb_ecmwf ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! ECMWF 453 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 454 CASE DEFAULT 455 CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 456 END SELECT 457 458 ! ! Compute true air density : 459 IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN ! At zu: (probably useless to remove zrho*grav*rn_zu from SLP...) 460 zrhoa(:,:) = rho_air( t_zu(:,:) , q_zu(:,:) , sf(jp_slp)%fnow(:,:,1) ) 461 ELSE ! At zt: 462 zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 463 END IF 464 465 !! CALL iom_put( "Cd_oce", Cd_atm) ! output value of pure ocean-atm. transfer coef. 466 !! CALL iom_put( "Ch_oce", Ch_atm) ! output value of pure ocean-atm. transfer coef. 467 468 DO jj = 1, jpj ! tau module, i and j component 469 DO ji = 1, jpi 470 zztmp = zrhoa(ji,jj) * zU_zu(ji,jj) * Cd_atm(ji,jj) ! using bulk wind speed 471 taum (ji,jj) = zztmp * wndm (ji,jj) 472 zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 473 zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 474 END DO 475 END DO 476 477 ! ! add the HF tau contribution to the wind stress module 478 IF( lhftau ) taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1) 479 480 CALL iom_put( "taum_oce", taum ) ! output wind stress module 481 482 ! ... utau, vtau at U- and V_points, resp. 483 ! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines 484 ! Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves 485 DO jj = 1, jpjm1 486 DO ji = 1, fs_jpim1 487 utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj ) ) & 488 & * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1)) 489 vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji ,jj+1) ) & 490 & * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1)) 491 END DO 492 END DO 493 CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. ) 494 495 ! Turbulent fluxes over ocean 496 ! ----------------------------- 497 498 ! zqla used as temporary array, for rho*U (common term of bulk formulae): 499 zqla(:,:) = zrhoa(:,:) * zU_zu(:,:) * tmask(:,:,1) 500 501 IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN 502 !! q_air and t_air are given at 10m (wind reference height) 503 zevap(:,:) = rn_efac*MAX( 0._wp, zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed 504 zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:) ) ! Sensible Heat, using bulk wind speed 505 ELSE 506 !! q_air and t_air are not given at 10m (wind reference height) 507 ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!! 508 zevap(:,:) = rn_efac*MAX( 0._wp, zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - q_zu(:,:) ) ) ! Evaporation, using bulk wind speed 509 zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) ) ! Sensible Heat, using bulk wind speed 510 ENDIF 511 512 zqla(:,:) = L_vap(zst(:,:)) * zevap(:,:) ! Latent Heat flux 513 514 515 IF(ln_ctl) THEN 516 CALL prt_ctl( tab2d_1=zqla , clinfo1=' blk_oce: zqla : ', tab2d_2=Ce_atm , clinfo2=' Ce_oce : ' ) 517 CALL prt_ctl( tab2d_1=zqsb , clinfo1=' blk_oce: zqsb : ', tab2d_2=Ch_atm , clinfo2=' Ch_oce : ' ) 518 CALL prt_ctl( tab2d_1=zqlw , clinfo1=' blk_oce: zqlw : ', tab2d_2=qsr, clinfo2=' qsr : ' ) 519 CALL prt_ctl( tab2d_1=zsq , clinfo1=' blk_oce: zsq : ', tab2d_2=zst, clinfo2=' zst : ' ) 520 CALL prt_ctl( tab2d_1=utau , clinfo1=' blk_oce: utau : ', mask1=umask, & 521 & tab2d_2=vtau , clinfo2= ' vtau : ', mask2=vmask ) 522 CALL prt_ctl( tab2d_1=wndm , clinfo1=' blk_oce: wndm : ') 523 CALL prt_ctl( tab2d_1=zst , clinfo1=' blk_oce: zst : ') 524 ENDIF 525 526 ! ----------------------------------------------------------------------------- ! 527 ! III Total FLUXES ! 528 ! ----------------------------------------------------------------------------- ! 529 ! 530 emp (:,:) = ( zevap(:,:) & ! mass flux (evap. - precip.) 531 & - sf(jp_prec)%fnow(:,:,1) * rn_pfac ) * tmask(:,:,1) 532 ! 533 qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:) & ! Downward Non Solar 534 & - sf(jp_snow)%fnow(:,:,1) * rn_pfac * rLfus & ! remove latent melting heat for solid precip 535 & - zevap(:,:) * pst(:,:) * rcp & ! remove evap heat content at SST 536 & + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) ) * rn_pfac & ! add liquid precip heat content at Tair 537 & * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp & 538 & + sf(jp_snow)%fnow(:,:,1) * rn_pfac & ! add solid precip heat content at min(Tair,Tsnow) 539 & * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi 849 ! 850 emp (:,:) = ( pevp(:,:) & ! mass flux (evap. - precip.) 851 & - pprec(:,:) * rn_pfac ) * tmask(:,:,1) 852 ! 853 qns(:,:) = zqlw(:,:) + psen(:,:) + zqla(:,:) & ! Downward Non Solar 854 & - psnow(:,:) * rn_pfac * rLfus & ! remove latent melting heat for solid precip 855 & - pevp(:,:) * ptsk(:,:) * rcp & ! remove evap heat content at SST 856 & + ( pprec(:,:) - psnow(:,:) ) * rn_pfac & ! add liquid precip heat content at Tair 857 & * ( ptair(:,:) - rt0 ) * rcp & 858 & + psnow(:,:) * rn_pfac & ! add solid precip heat content at min(Tair,Tsnow) 859 & * ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi 540 860 qns(:,:) = qns(:,:) * tmask(:,:,1) 541 861 ! 542 862 #if defined key_si3 543 qns_oce(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)! non solar without emp (only needed by SI3)863 qns_oce(:,:) = zqlw(:,:) + psen(:,:) + zqla(:,:) ! non solar without emp (only needed by SI3) 544 864 qsr_oce(:,:) = qsr(:,:) 545 865 #endif 546 866 ! 867 CALL iom_put( "rho_air" , rhoa*tmask(:,:,1) ) ! output air density [kg/m^3] 868 CALL iom_put( "evap_oce" , pevp ) ! evaporation 869 CALL iom_put( "qlw_oce" , zqlw ) ! output downward longwave heat over the ocean 870 CALL iom_put( "qsb_oce" , psen ) ! output downward sensible heat over the ocean 871 CALL iom_put( "qla_oce" , zqla ) ! output downward latent heat over the ocean 872 tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s] 873 sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s] 874 CALL iom_put( 'snowpre', sprecip ) ! Snow 875 CALL iom_put( 'precip' , tprecip ) ! Total precipitation 876 ! 547 877 IF ( nn_ice == 0 ) THEN 548 CALL iom_put( "qlw_oce" , zqlw ) ! output downward longwave heat over the ocean 549 CALL iom_put( "qsb_oce" , - zqsb ) ! output downward sensible heat over the ocean 550 CALL iom_put( "qla_oce" , - zqla ) ! output downward latent heat over the ocean 551 CALL iom_put( "qemp_oce", qns-zqlw+zqsb+zqla ) ! output downward heat content of E-P over the ocean 552 CALL iom_put( "qns_oce" , qns ) ! output downward non solar heat over the ocean 553 CALL iom_put( "qsr_oce" , qsr ) ! output downward solar heat over the ocean 554 CALL iom_put( "qt_oce" , qns+qsr ) ! output total downward heat over the ocean 555 tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s] 556 sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s] 557 CALL iom_put( 'snowpre', sprecip ) ! Snow 558 CALL iom_put( 'precip' , tprecip ) ! Total precipitation 559 ENDIF 560 ! 561 IF(ln_ctl) THEN 562 CALL prt_ctl(tab2d_1=zqsb , clinfo1=' blk_oce: zqsb : ', tab2d_2=zqlw , clinfo2=' zqlw : ') 563 CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce: zqla : ', tab2d_2=qsr , clinfo2=' qsr : ') 564 CALL prt_ctl(tab2d_1=pst , clinfo1=' blk_oce: pst : ', tab2d_2=emp , clinfo2=' emp : ') 565 CALL prt_ctl(tab2d_1=utau , clinfo1=' blk_oce: utau : ', mask1=umask, & 566 & tab2d_2=vtau , clinfo2= ' vtau : ' , mask2=vmask ) 567 ENDIF 568 ! 569 END SUBROUTINE blk_oce 570 571 572 573 FUNCTION rho_air( ptak, pqa, pslp ) 574 !!------------------------------------------------------------------------------- 575 !! *** FUNCTION rho_air *** 576 !! 577 !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere 578 !! 579 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 580 !!------------------------------------------------------------------------------- 581 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K] 582 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg] 583 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp ! pressure in [Pa] 584 REAL(wp), DIMENSION(jpi,jpj) :: rho_air ! density of moist air [kg/m^3] 585 !!------------------------------------------------------------------------------- 586 ! 587 rho_air = pslp / ( R_dry*ptak * ( 1._wp + rctv0*pqa ) ) 588 ! 589 END FUNCTION rho_air 590 591 592 FUNCTION cp_air( pqa ) 593 !!------------------------------------------------------------------------------- 594 !! *** FUNCTION cp_air *** 595 !! 596 !! ** Purpose : Compute specific heat (Cp) of moist air 597 !! 598 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 599 !!------------------------------------------------------------------------------- 600 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg] 601 REAL(wp), DIMENSION(jpi,jpj) :: cp_air ! specific heat of moist air [J/K/kg] 602 !!------------------------------------------------------------------------------- 603 ! 604 Cp_air = Cp_dry + Cp_vap * pqa 605 ! 606 END FUNCTION cp_air 607 608 609 FUNCTION q_sat( ptak, pslp ) 610 !!---------------------------------------------------------------------------------- 611 !! *** FUNCTION q_sat *** 612 !! 613 !! ** Purpose : Specific humidity at saturation in [kg/kg] 614 !! Based on accurate estimate of "e_sat" 615 !! aka saturation water vapor (Goff, 1957) 616 !! 617 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 618 !!---------------------------------------------------------------------------------- 619 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K] 620 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp ! sea level atmospheric pressure [Pa] 621 REAL(wp), DIMENSION(jpi,jpj) :: q_sat ! Specific humidity at saturation [kg/kg] 622 ! 623 INTEGER :: ji, jj ! dummy loop indices 624 REAL(wp) :: ze_sat, ztmp ! local scalar 625 !!---------------------------------------------------------------------------------- 626 ! 627 DO jj = 1, jpj 628 DO ji = 1, jpi 629 ! 630 ztmp = rt0 / ptak(ji,jj) 631 ! 632 ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957) 633 ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0) & 634 & + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) ) & 635 & + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614 ) 636 ! 637 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] 638 ! 639 END DO 640 END DO 641 ! 642 END FUNCTION q_sat 643 644 645 FUNCTION gamma_moist( ptak, pqa ) 646 !!---------------------------------------------------------------------------------- 647 !! *** FUNCTION gamma_moist *** 648 !! 649 !! ** Purpose : Compute the moist adiabatic lapse-rate. 650 !! => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate 651 !! => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html 652 !! 653 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 654 !!---------------------------------------------------------------------------------- 655 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K] 656 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! specific humidity [kg/kg] 657 REAL(wp), DIMENSION(jpi,jpj) :: gamma_moist ! moist adiabatic lapse-rate 658 ! 659 INTEGER :: ji, jj ! dummy loop indices 660 REAL(wp) :: zrv, ziRT ! local scalar 661 !!---------------------------------------------------------------------------------- 662 ! 663 DO jj = 1, jpj 664 DO ji = 1, jpi 665 zrv = pqa(ji,jj) / (1. - pqa(ji,jj)) 666 ziRT = 1. / (R_dry*ptak(ji,jj)) ! 1/RT 667 gamma_moist(ji,jj) = grav * ( 1. + rLevap*zrv*ziRT ) / ( Cp_dry + rLevap*rLevap*zrv*reps0*ziRT/ptak(ji,jj) ) 668 END DO 669 END DO 670 ! 671 END FUNCTION gamma_moist 672 673 674 FUNCTION L_vap( psst ) 675 !!--------------------------------------------------------------------------------- 676 !! *** FUNCTION L_vap *** 677 !! 678 !! ** Purpose : Compute the latent heat of vaporization of water from temperature 679 !! 680 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 681 !!---------------------------------------------------------------------------------- 682 REAL(wp), DIMENSION(jpi,jpj) :: L_vap ! latent heat of vaporization [J/kg] 683 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst ! water temperature [K] 684 !!---------------------------------------------------------------------------------- 685 ! 686 L_vap = ( 2.501 - 0.00237 * ( psst(:,:) - rt0) ) * 1.e6 687 ! 688 END FUNCTION L_vap 878 CALL iom_put( "qemp_oce" , qns-zqlw-psen-zqla ) ! output downward heat content of E-P over the ocean 879 CALL iom_put( "qns_oce" , qns ) ! output downward non solar heat over the ocean 880 CALL iom_put( "qsr_oce" , qsr ) ! output downward solar heat over the ocean 881 CALL iom_put( "qt_oce" , qns+qsr ) ! output total downward heat over the ocean 882 ENDIF 883 ! 884 IF(sn_cfctl%l_prtctl) THEN 885 CALL prt_ctl(tab2d_1=zqlw , clinfo1=' blk_oce_2: zqlw : ') 886 CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce_2: zqla : ', tab2d_2=qsr , clinfo2=' qsr : ') 887 CALL prt_ctl(tab2d_1=emp , clinfo1=' blk_oce_2: emp : ') 888 ENDIF 889 ! 890 END SUBROUTINE blk_oce_2 891 689 892 690 893 SUBROUTINE iom_blk_wgt_init( cdname, ld_tmppatch ) … … 857 1060 !! 'key_si3' SI3 sea-ice model 858 1061 !!---------------------------------------------------------------------- 859 !! blk_ice_ tau: provide the air-ice stress860 !! blk_ice_ flx: provide the heat and mass fluxes at air-ice interface1062 !! blk_ice_1 : provide the air-ice stress 1063 !! blk_ice_2 : provide the heat and mass fluxes at air-ice interface 861 1064 !! blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux) 862 1065 !! Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag 863 !! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 1066 !! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 864 1067 !!---------------------------------------------------------------------- 865 1068 866 SUBROUTINE blk_ice_tau 867 !!--------------------------------------------------------------------- 868 !! *** ROUTINE blk_ice_tau *** 1069 SUBROUTINE blk_ice_1( pwndi, pwndj, ptair, phumi, pslp , puice, pvice, ptsui, & ! inputs 1070 & putaui, pvtaui, pseni, pevpi, pssqi, pcd_dui ) ! optional outputs 1071 !!--------------------------------------------------------------------- 1072 !! *** ROUTINE blk_ice_1 *** 869 1073 !! 870 1074 !! ** Purpose : provide the surface boundary condition over sea-ice … … 874 1078 !! NB: ice drag coefficient is assumed to be a constant 875 1079 !!--------------------------------------------------------------------- 1080 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pslp ! sea-level pressure [Pa] 1081 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndi ! atmospheric wind at T-point [m/s] 1082 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndj ! atmospheric wind at T-point [m/s] 1083 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: ptair ! atmospheric wind at T-point [m/s] 1084 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: phumi ! atmospheric wind at T-point [m/s] 1085 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: puice ! sea-ice velocity on I or C grid [m/s] 1086 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pvice ! " 1087 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: ptsui ! sea-ice surface temperature [K] 1088 REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: putaui ! if ln_blk 1089 REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pvtaui ! if ln_blk 1090 REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pseni ! if ln_abl 1091 REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pevpi ! if ln_abl 1092 REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pssqi ! if ln_abl 1093 REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pcd_dui ! if ln_abl 1094 ! 876 1095 INTEGER :: ji, jj ! dummy loop indices 877 REAL(wp) :: zwndi_f , zwndj_f, zwnorm_f ! relative wind module and components at F-point 878 REAL(wp) :: zwndi_t , zwndj_t ! relative wind components at T-point 879 REAL(wp), DIMENSION(jpi,jpj) :: zrhoa ! transfer coefficient for momentum (tau) 880 !!--------------------------------------------------------------------- 881 ! 882 ! set transfer coefficients to default sea-ice values 883 Cd_atm(:,:) = Cd_ice 884 Ch_atm(:,:) = Cd_ice 885 Ce_atm(:,:) = Cd_ice 886 887 wndm_ice(:,:) = 0._wp !!gm brutal.... 1096 REAL(wp) :: zootm_su ! sea-ice surface mean temperature 1097 REAL(wp) :: zztmp1, zztmp2 ! temporary arrays 1098 REAL(wp), DIMENSION(jpi,jpj) :: zcd_dui ! transfer coefficient for momentum (tau) 1099 !!--------------------------------------------------------------------- 1100 ! 888 1101 889 1102 ! ------------------------------------------------------------ ! … … 891 1104 ! ------------------------------------------------------------ ! 892 1105 ! C-grid ice dynamics : U & V-points (same as ocean) 893 DO jj = 2, jpjm1 894 DO ji = fs_2, fs_jpim1 ! vect. opt. 895 zwndi_t = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) ) 896 zwndj_t = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji ,jj-1) + v_ice(ji,jj) ) ) 897 wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 898 END DO 899 END DO 900 CALL lbc_lnk( 'sbcblk', wndm_ice, 'T', 1. ) 1106 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 1107 wndm_ice(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) ) 1108 END_2D 901 1109 ! 902 1110 ! Make ice-atm. drag dependent on ice concentration 903 1111 IF ( ln_Cd_L12 ) THEN ! calculate new drag from Lupkes(2012) equations 904 CALL Cdn10_Lupkes2012( Cd_atm ) 905 Ch_atm(:,:) = Cd_atm(:,:) ! momentum and heat transfer coef. are considered identical 1112 CALL Cdn10_Lupkes2012( Cd_ice ) 1113 Ch_ice(:,:) = Cd_ice(:,:) ! momentum and heat transfer coef. are considered identical 1114 Ce_ice(:,:) = Cd_ice(:,:) 906 1115 ELSEIF( ln_Cd_L15 ) THEN ! calculate new drag from Lupkes(2015) equations 907 CALL Cdn10_Lupkes2015( Cd_atm, Ch_atm ) 908 ENDIF 909 910 !! CALL iom_put( "Cd_ice", Cd_atm) ! output value of pure ice-atm. transfer coef. 911 !! CALL iom_put( "Ch_ice", Ch_atm) ! output value of pure ice-atm. transfer coef. 912 1116 CALL Cdn10_Lupkes2015( ptsui, pslp, Cd_ice, Ch_ice ) 1117 Ce_ice(:,:) = Ch_ice(:,:) ! sensible and latent heat transfer coef. are considered identical 1118 ENDIF 1119 1120 IF( iom_use('Cd_ice') ) CALL iom_put("Cd_ice", Cd_ice) 1121 IF( iom_use('Ce_ice') ) CALL iom_put("Ce_ice", Ce_ice) 1122 IF( iom_use('Ch_ice') ) CALL iom_put("Ch_ice", Ch_ice) 1123 913 1124 ! local scalars ( place there for vector optimisation purposes) 914 ! Computing density of air! Way denser that 1.2 over sea-ice !!! 915 zrhoa (:,:) = rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) 916 917 !!gm brutal.... 918 utau_ice (:,:) = 0._wp 919 vtau_ice (:,:) = 0._wp 920 !!gm end 921 922 ! ------------------------------------------------------------ ! 923 ! Wind stress relative to the moving ice ( U10m - U_ice ) ! 924 ! ------------------------------------------------------------ ! 925 ! C-grid ice dynamics : U & V-points (same as ocean) 926 DO jj = 2, jpjm1 927 DO ji = fs_2, fs_jpim1 ! vect. opt. 928 utau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji+1,jj ) + wndm_ice(ji,jj) ) & 929 & * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) ) 930 vtau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji,jj+1 ) + wndm_ice(ji,jj) ) & 931 & * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) ) 932 END DO 933 END DO 934 CALL lbc_lnk_multi( 'sbcblk', utau_ice, 'U', -1., vtau_ice, 'V', -1. ) 935 ! 936 ! 937 IF(ln_ctl) THEN 938 CALL prt_ctl(tab2d_1=utau_ice , clinfo1=' blk_ice: utau_ice : ', tab2d_2=vtau_ice , clinfo2=' vtau_ice : ') 939 CALL prt_ctl(tab2d_1=wndm_ice , clinfo1=' blk_ice: wndm_ice : ') 940 ENDIF 941 ! 942 END SUBROUTINE blk_ice_tau 943 944 945 SUBROUTINE blk_ice_flx( ptsu, phs, phi, palb ) 946 !!--------------------------------------------------------------------- 947 !! *** ROUTINE blk_ice_flx *** 1125 zcd_dui(:,:) = wndm_ice(:,:) * Cd_ice(:,:) 1126 1127 IF( ln_blk ) THEN 1128 ! ---------------------------------------------------- ! 1129 ! Wind stress relative to nonmoving ice ( U10m ) ! 1130 ! ---------------------------------------------------- ! 1131 ! supress moving ice in wind stress computation as we don't know how to do it properly... 1132 DO_2D( 0, 1, 0, 1 ) ! at T point 1133 putaui(ji,jj) = rhoa(ji,jj) * zcd_dui(ji,jj) * pwndi(ji,jj) 1134 pvtaui(ji,jj) = rhoa(ji,jj) * zcd_dui(ji,jj) * pwndj(ji,jj) 1135 END_2D 1136 ! 1137 DO_2D( 0, 0, 0, 0 ) ! U & V-points (same as ocean). 1138 ! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and rheology 1139 zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj ,1) ) 1140 zztmp2 = 0.5_wp * ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji ,jj+1,1) ) 1141 putaui(ji,jj) = zztmp1 * ( putaui(ji,jj) + putaui(ji+1,jj ) ) 1142 pvtaui(ji,jj) = zztmp2 * ( pvtaui(ji,jj) + pvtaui(ji ,jj+1) ) 1143 END_2D 1144 CALL lbc_lnk_multi( 'sbcblk', putaui, 'U', -1._wp, pvtaui, 'V', -1._wp ) 1145 ! 1146 IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=putaui , clinfo1=' blk_ice: putaui : ' & 1147 & , tab2d_2=pvtaui , clinfo2=' pvtaui : ' ) 1148 ELSE ! ln_abl 1149 zztmp1 = 11637800.0_wp 1150 zztmp2 = -5897.8_wp 1151 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 1152 pcd_dui(ji,jj) = zcd_dui (ji,jj) 1153 pseni (ji,jj) = wndm_ice(ji,jj) * Ch_ice(ji,jj) 1154 pevpi (ji,jj) = wndm_ice(ji,jj) * Ce_ice(ji,jj) 1155 zootm_su = zztmp2 / ptsui(ji,jj) ! ptsui is in K (it can't be zero ??) 1156 pssqi (ji,jj) = zztmp1 * EXP( zootm_su ) / rhoa(ji,jj) 1157 END_2D 1158 ENDIF 1159 ! 1160 IF(sn_cfctl%l_prtctl) CALL prt_ctl(tab2d_1=wndm_ice , clinfo1=' blk_ice: wndm_ice : ') 1161 ! 1162 END SUBROUTINE blk_ice_1 1163 1164 1165 SUBROUTINE blk_ice_2( ptsu, phs, phi, palb, ptair, phumi, pslp, pqlw, pprec, psnow ) 1166 !!--------------------------------------------------------------------- 1167 !! *** ROUTINE blk_ice_2 *** 948 1168 !! 949 1169 !! ** Purpose : provide the heat and mass fluxes at air-ice interface … … 955 1175 !! caution : the net upward water flux has with mm/day unit 956 1176 !!--------------------------------------------------------------------- 957 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: ptsu ! sea ice surface temperature 1177 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: ptsu ! sea ice surface temperature [K] 958 1178 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phs ! snow thickness 959 1179 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phi ! ice thickness 960 1180 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: palb ! ice albedo (all skies) 1181 REAL(wp), DIMENSION(:,: ), INTENT(in) :: ptair 1182 REAL(wp), DIMENSION(:,: ), INTENT(in) :: phumi 1183 REAL(wp), DIMENSION(:,: ), INTENT(in) :: pslp 1184 REAL(wp), DIMENSION(:,: ), INTENT(in) :: pqlw 1185 REAL(wp), DIMENSION(:,: ), INTENT(in) :: pprec 1186 REAL(wp), DIMENSION(:,: ), INTENT(in) :: psnow 961 1187 !! 962 1188 INTEGER :: ji, jj, jl ! dummy loop indices 963 1189 REAL(wp) :: zst3 ! local variable 964 1190 REAL(wp) :: zcoef_dqlw, zcoef_dqla ! - - 965 REAL(wp) :: zztmp, z 1_rLsub! - -1191 REAL(wp) :: zztmp, zztmp2, z1_rLsub ! - - 966 1192 REAL(wp) :: zfr1, zfr2 ! local variables 967 1193 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_st ! inverse of surface temperature … … 971 1197 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqsb ! sensible heat sensitivity over ice 972 1198 REAL(wp), DIMENSION(jpi,jpj) :: zevap, zsnw ! evaporation and snw distribution after wind blowing (SI3) 973 REAL(wp), DIMENSION(jpi,jpj) :: zrhoa 974 !!--------------------------------------------------------------------- 975 ! 976 zcoef_dqlw = 4.0 * 0.95 * Stef ! local scalars 977 zcoef_dqla = -Ls * 11637800. * (-5897.8) 978 ! 979 zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 1199 REAL(wp), DIMENSION(jpi,jpj) :: zqair ! specific humidity of air at z=rn_zqt [kg/kg] !LB 1200 REAL(wp), DIMENSION(jpi,jpj) :: ztmp, ztmp2 1201 !!--------------------------------------------------------------------- 1202 ! 1203 zcoef_dqlw = 4._wp * 0.95_wp * stefan ! local scalars 1204 zcoef_dqla = -rLsub * 11637800._wp * (-5897.8_wp) !LB: BAD! 1205 ! 1206 SELECT CASE( nhumi ) 1207 CASE( np_humi_sph ) 1208 zqair(:,:) = phumi(:,:) ! what we read in file is already a spec. humidity! 1209 CASE( np_humi_dpt ) 1210 zqair(:,:) = q_sat( phumi(:,:), pslp ) 1211 CASE( np_humi_rlh ) 1212 zqair(:,:) = q_air_rh( 0.01_wp*phumi(:,:), ptair(:,:), pslp(:,:) ) !LB: 0.01 => RH is % percent in file 1213 END SELECT 980 1214 ! 981 1215 zztmp = 1. / ( 1. - albo ) 982 WHERE( ptsu(:,:,:) /= 0._wp ) ; z1_st(:,:,:) = 1._wp / ptsu(:,:,:) 983 ELSEWHERE ; z1_st(:,:,:) = 0._wp 1216 WHERE( ptsu(:,:,:) /= 0._wp ) 1217 z1_st(:,:,:) = 1._wp / ptsu(:,:,:) 1218 ELSEWHERE 1219 z1_st(:,:,:) = 0._wp 984 1220 END WHERE 985 1221 ! ! ========================== ! … … 995 1231 qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj) 996 1232 ! Long Wave (lw) 997 z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef* ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)1233 z_qlw(ji,jj,jl) = 0.95 * ( pqlw(ji,jj) - stefan * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 998 1234 ! lw sensitivity 999 1235 z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3 … … 1003 1239 ! ----------------------------! 1004 1240 1005 ! ... turbulent heat fluxes with Ch_ atm recalculated in blk_ice_tau1241 ! ... turbulent heat fluxes with Ch_ice recalculated in blk_ice_1 1006 1242 ! Sensible Heat 1007 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))1243 z_qsb(ji,jj,jl) = rhoa(ji,jj) * rCp_air * Ch_ice(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - ptair(ji,jj)) 1008 1244 ! Latent Heat 1009 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls * Ch_atm(ji,jj) * wndm_ice(ji,jj) * & 1010 & ( 11637800. * EXP( -5897.8 * z1_st(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) ) 1245 zztmp2 = EXP( -5897.8 * z1_st(ji,jj,jl) ) 1246 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa(ji,jj) * rLsub * Ce_ice(ji,jj) * wndm_ice(ji,jj) * & 1247 & ( 11637800. * zztmp2 / rhoa(ji,jj) - zqair(ji,jj) ) ) 1011 1248 ! Latent heat sensitivity for ice (Dqla/Dt) 1012 1249 IF( qla_ice(ji,jj,jl) > 0._wp ) THEN 1013 dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * C h_atm(ji,jj) * wndm_ice(ji,jj) * &1014 & z1_st(ji,jj,jl) *z1_st(ji,jj,jl) * EXP(-5897.8 * z1_st(ji,jj,jl))1250 dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ce_ice(ji,jj) * wndm_ice(ji,jj) * & 1251 & z1_st(ji,jj,jl) * z1_st(ji,jj,jl) * zztmp2 1015 1252 ELSE 1016 1253 dqla_ice(ji,jj,jl) = 0._wp … … 1018 1255 1019 1256 ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice) 1020 z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj)1257 z_dqsb(ji,jj,jl) = rhoa(ji,jj) * rCp_air * Ch_ice(ji,jj) * wndm_ice(ji,jj) 1021 1258 1022 1259 ! ----------------------------! … … 1033 1270 END DO 1034 1271 ! 1035 tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! total precipitation [kg/m2/s]1036 sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! solid precipitation [kg/m2/s]1037 CALL iom_put( 'snowpre', sprecip ) 1038 CALL iom_put( 'precip' , tprecip ) 1272 tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1) ! total precipitation [kg/m2/s] 1273 sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1) ! solid precipitation [kg/m2/s] 1274 CALL iom_put( 'snowpre', sprecip ) ! Snow precipitation 1275 CALL iom_put( 'precip' , tprecip ) ! Total precipitation 1039 1276 1040 1277 ! --- evaporation --- ! … … 1042 1279 evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_rLsub ! sublimation 1043 1280 devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_rLsub ! d(sublimation)/dT 1044 zevap (:,:) = rn_efac * ( emp(:,:) + tprecip(:,:) ) ! evaporation over ocean1281 zevap (:,:) = emp(:,:) + tprecip(:,:) ! evaporation over ocean !LB: removed rn_efac here, correct??? 1045 1282 1046 1283 ! --- evaporation minus precipitation --- ! … … 1053 1290 ! --- heat flux associated with emp --- ! 1054 1291 qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp & ! evap at sst 1055 & + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp& ! liquid precip at Tair1292 & + ( tprecip(:,:) - sprecip(:,:) ) * ( ptair(:,:) - rt0 ) * rcp & ! liquid precip at Tair 1056 1293 & + sprecip(:,:) * ( 1._wp - zsnw ) * & ! solid precip at min(Tair,Tsnow) 1057 & ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )1294 & ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 1058 1295 qemp_ice(:,:) = sprecip(:,:) * zsnw * & ! solid precip (only) 1059 & ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )1296 & ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 1060 1297 1061 1298 ! --- total solar and non solar fluxes --- ! … … 1065 1302 1066 1303 ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- ! 1067 qprec_ice(:,:) = rhos * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )1304 qprec_ice(:,:) = rhos * ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 1068 1305 1069 1306 ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- 1070 1307 DO jl = 1, jpl 1071 1308 qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * rcpi * tmask(:,:,1) ) 1072 ! ! But we do not have Tice => consider it at 0degC => evap=0 1309 ! ! But we do not have Tice => consider it at 0degC => evap=0 1073 1310 END DO 1074 1311 … … 1077 1314 zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) ! zfr2 such that zfr1 + zfr2 to equal 1 1078 1315 ! 1079 WHERE ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) < 0.1_wp ) ! linear decrease from hi=0 to 10cm 1316 WHERE ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) < 0.1_wp ) ! linear decrease from hi=0 to 10cm 1080 1317 qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) ) 1081 1318 ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp ) ! constant (zfr1) when hi>10cm 1082 1319 qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1 1083 1320 ELSEWHERE ! zero when hs>0 1084 qtr_ice_top(:,:,:) = 0._wp 1321 qtr_ice_top(:,:,:) = 0._wp 1085 1322 END WHERE 1086 1323 ! 1087 IF(ln_ctl) THEN 1324 1325 IF( iom_use('evap_ao_cea') .OR. iom_use('hflx_evap_cea') ) THEN 1326 ztmp(:,:) = zevap(:,:) * ( 1._wp - at_i_b(:,:) ) 1327 IF( iom_use('evap_ao_cea' ) ) CALL iom_put( 'evap_ao_cea' , ztmp(:,:) * tmask(:,:,1) ) ! ice-free oce evap (cell average) 1328 IF( iom_use('hflx_evap_cea') ) CALL iom_put( 'hflx_evap_cea', ztmp(:,:) * sst_m(:,:) * rcp * tmask(:,:,1) ) ! heat flux from evap (cell average) 1329 ENDIF 1330 IF( iom_use('hflx_rain_cea') ) THEN 1331 ztmp(:,:) = rcp * ( SUM( (ptsu-rt0) * a_i_b, dim=3 ) + sst_m(:,:) * ( 1._wp - at_i_b(:,:) ) ) 1332 IF( iom_use('hflx_rain_cea') ) CALL iom_put( 'hflx_rain_cea', ( tprecip(:,:) - sprecip(:,:) ) * ztmp(:,:) ) ! heat flux from rain (cell average) 1333 ENDIF 1334 IF( iom_use('hflx_snow_cea') .OR. iom_use('hflx_snow_ao_cea') .OR. iom_use('hflx_snow_ai_cea') ) THEN 1335 WHERE( SUM( a_i_b, dim=3 ) > 1.e-10 ) 1336 ztmp(:,:) = rcpi * SUM( (ptsu-rt0) * a_i_b, dim=3 ) / SUM( a_i_b, dim=3 ) 1337 ELSEWHERE 1338 ztmp(:,:) = rcp * sst_m(:,:) 1339 ENDWHERE 1340 ztmp2(:,:) = sprecip(:,:) * ( ztmp(:,:) - rLfus ) 1341 IF( iom_use('hflx_snow_cea') ) CALL iom_put('hflx_snow_cea' , ztmp2(:,:) ) ! heat flux from snow (cell average) 1342 IF( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea', ztmp2(:,:) * ( 1._wp - zsnw(:,:) ) ) ! heat flux from snow (over ocean) 1343 IF( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea', ztmp2(:,:) * zsnw(:,:) ) ! heat flux from snow (over ice) 1344 ENDIF 1345 ! 1346 IF(sn_cfctl%l_prtctl) THEN 1088 1347 CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice : ', tab3d_2=z_qsb , clinfo2=' z_qsb : ', kdim=jpl) 1089 1348 CALL prt_ctl(tab3d_1=z_qlw , clinfo1=' blk_ice: z_qlw : ', tab3d_2=dqla_ice, clinfo2=' dqla_ice : ', kdim=jpl) … … 1094 1353 ENDIF 1095 1354 ! 1096 END SUBROUTINE blk_ice_ flx1097 1355 END SUBROUTINE blk_ice_2 1356 1098 1357 1099 1358 SUBROUTINE blk_ice_qcn( ld_virtual_itd, ptsu, ptb, phs, phi ) … … 1104 1363 !! to force sea ice / snow thermodynamics 1105 1364 !! in the case conduction flux is emulated 1106 !! 1365 !! 1107 1366 !! ** Method : compute surface energy balance assuming neglecting heat storage 1108 1367 !! following the 0-layer Semtner (1976) approach … … 1129 1388 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zgfac ! enhanced conduction factor 1130 1389 !!--------------------------------------------------------------------- 1131 1390 1132 1391 ! -------------------------------------! 1133 1392 ! I Enhanced conduction factor ! … … 1137 1396 ! 1138 1397 zgfac(:,:,:) = 1._wp 1139 1398 1140 1399 IF( ld_virtual_itd ) THEN 1141 1400 ! … … 1143 1402 zfac2 = EXP(1._wp) * 0.5_wp * zepsilon 1144 1403 zfac3 = 2._wp / zepsilon 1145 ! 1146 DO jl = 1, jpl 1147 DO jj = 1 , jpj 1148 DO ji = 1, jpi 1149 zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac ! Effective thickness 1150 IF( zhe >= zfac2 ) zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor 1151 END DO 1152 END DO 1404 ! 1405 DO jl = 1, jpl 1406 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 1407 zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac ! Effective thickness 1408 IF( zhe >= zfac2 ) zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor 1409 END_2D 1153 1410 END DO 1154 ! 1155 ENDIF 1156 1411 ! 1412 ENDIF 1413 1157 1414 ! -------------------------------------------------------------! 1158 1415 ! II Surface temperature and conduction flux ! … … 1162 1419 ! 1163 1420 DO jl = 1, jpl 1164 DO jj = 1 , jpj 1165 DO ji = 1, jpi 1166 ! 1167 zkeff_h = zfac * zgfac(ji,jj,jl) / & ! Effective conductivity of the snow-ice system divided by thickness 1168 & ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) ) 1169 ztsu = ptsu(ji,jj,jl) ! Store current iteration temperature 1170 ztsu0 = ptsu(ji,jj,jl) ! Store initial surface temperature 1171 zqa0 = qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) ! Net initial atmospheric heat flux 1172 ! 1173 DO iter = 1, nit ! --- Iterative loop 1174 zqc = zkeff_h * ( ztsu - ptb(ji,jj) ) ! Conduction heat flux through snow-ice system (>0 downwards) 1175 zqnet = zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc ! Surface energy budget 1176 ztsu = ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h ) ! Temperature update 1177 END DO 1178 ! 1179 ptsu (ji,jj,jl) = MIN( rt0, ztsu ) 1180 qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) ) 1181 qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) 1182 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) ) & 1183 & * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) ) 1184 1185 ! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- ! 1186 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) 1187 1421 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 1422 ! 1423 zkeff_h = zfac * zgfac(ji,jj,jl) / & ! Effective conductivity of the snow-ice system divided by thickness 1424 & ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) ) 1425 ztsu = ptsu(ji,jj,jl) ! Store current iteration temperature 1426 ztsu0 = ptsu(ji,jj,jl) ! Store initial surface temperature 1427 zqa0 = qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) ! Net initial atmospheric heat flux 1428 ! 1429 DO iter = 1, nit ! --- Iterative loop 1430 zqc = zkeff_h * ( ztsu - ptb(ji,jj) ) ! Conduction heat flux through snow-ice system (>0 downwards) 1431 zqnet = zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc ! Surface energy budget 1432 ztsu = ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h ) ! Temperature update 1188 1433 END DO 1189 END DO 1434 ! 1435 ptsu (ji,jj,jl) = MIN( rt0, ztsu ) 1436 qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) ) 1437 qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) 1438 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) ) & 1439 & * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) ) 1440 1441 ! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- ! 1442 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) 1443 1444 END_2D 1190 1445 ! 1191 END DO 1192 ! 1446 END DO 1447 ! 1193 1448 END SUBROUTINE blk_ice_qcn 1194 1195 1196 SUBROUTINE Cdn10_Lupkes2012( Cd )1449 1450 1451 SUBROUTINE Cdn10_Lupkes2012( pcd ) 1197 1452 !!---------------------------------------------------------------------- 1198 1453 !! *** ROUTINE Cdn10_Lupkes2012 *** 1199 1454 !! 1200 !! ** Purpose : Recompute the neutral air-ice drag referenced at 10m 1455 !! ** Purpose : Recompute the neutral air-ice drag referenced at 10m 1201 1456 !! to make it dependent on edges at leads, melt ponds and flows. 1202 1457 !! After some approximations, this can be resumed to a dependency 1203 1458 !! on ice concentration. 1204 !! 1459 !! 1205 1460 !! ** Method : The parameterization is taken from Lupkes et al. (2012) eq.(50) 1206 1461 !! with the highest level of approximation: level4, eq.(59) … … 1214 1469 !! 1215 1470 !! This new drag has a parabolic shape (as a function of A) starting at 1216 !! Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5 1471 !! Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5 1217 1472 !! and going down to Cdi(say 1.4e-3) for A=1 1218 1473 !! … … 1224 1479 !! 1225 1480 !!---------------------------------------------------------------------- 1226 REAL(wp), DIMENSION(:,:), INTENT(inout) :: Cd1481 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pcd 1227 1482 REAL(wp), PARAMETER :: zCe = 2.23e-03_wp 1228 1483 REAL(wp), PARAMETER :: znu = 1._wp … … 1239 1494 1240 1495 ! ice-atm drag 1241 Cd(:,:) =Cd_ice + & ! pure ice drag1242 & zCe * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp) ! change due to sea-ice morphology1243 1496 pcd(:,:) = rCd_ice + & ! pure ice drag 1497 & zCe * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp) ! change due to sea-ice morphology 1498 1244 1499 END SUBROUTINE Cdn10_Lupkes2012 1245 1500 1246 1501 1247 SUBROUTINE Cdn10_Lupkes2015( Cd, Ch )1502 SUBROUTINE Cdn10_Lupkes2015( ptm_su, pslp, pcd, pch ) 1248 1503 !!---------------------------------------------------------------------- 1249 1504 !! *** ROUTINE Cdn10_Lupkes2015 *** 1250 1505 !! 1251 1506 !! ** pUrpose : Alternative turbulent transfert coefficients formulation 1252 !! between sea-ice and atmosphere with distinct momentum 1253 !! and heat coefficients depending on sea-ice concentration 1507 !! between sea-ice and atmosphere with distinct momentum 1508 !! and heat coefficients depending on sea-ice concentration 1254 1509 !! and atmospheric stability (no meltponds effect for now). 1255 !! 1510 !! 1256 1511 !! ** Method : The parameterization is adapted from Lupkes et al. (2015) 1257 1512 !! and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme, 1258 1513 !! it considers specific skin and form drags (Andreas et al. 2010) 1259 !! to compute neutral transfert coefficients for both heat and 1514 !! to compute neutral transfert coefficients for both heat and 1260 1515 !! momemtum fluxes. Atmospheric stability effect on transfert 1261 1516 !! coefficient is also taken into account following Louis (1979). … … 1266 1521 !!---------------------------------------------------------------------- 1267 1522 ! 1268 REAL(wp), DIMENSION(:,:), INTENT(inout) :: Cd 1269 REAL(wp), DIMENSION(:,:), INTENT(inout) :: Ch 1270 REAL(wp), DIMENSION(jpi,jpj) :: ztm_su, zst, zqo_sat, zqi_sat 1523 REAL(wp), DIMENSION(:,:), INTENT(in ) :: ptm_su ! sea-ice surface temperature [K] 1524 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pslp ! sea-level pressure [Pa] 1525 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pcd ! momentum transfert coefficient 1526 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pch ! heat transfert coefficient 1527 REAL(wp), DIMENSION(jpi,jpj) :: zst, zqo_sat, zqi_sat 1271 1528 ! 1272 1529 ! ECHAM6 constants … … 1296 1553 !!---------------------------------------------------------------------- 1297 1554 1298 ! mean temperature1299 WHERE( at_i_b(:,:) > 1.e-20 ) ; ztm_su(:,:) = SUM( t_su(:,:,:) * a_i_b(:,:,:) , dim=3 ) / at_i_b(:,:)1300 ELSEWHERE ; ztm_su(:,:) = rt01301 ENDWHERE1302 1303 1555 ! Momentum Neutral Transfert Coefficients (should be a constant) 1304 1556 zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2 ! Eq. 40 1305 1557 zCdn_skin_ice = ( vkarmn / LOG( rn_zu / z0_skin_ice + 1._wp ) )**2 ! Eq. 7 1306 zCdn_ice = zCdn_skin_ice ! Eq. 7 (cf Lupkes email for details)1558 zCdn_ice = zCdn_skin_ice ! Eq. 7 1307 1559 !zCdn_ice = 1.89e-3 ! old ECHAM5 value (cf Eq. 32) 1308 1560 1309 1561 ! Heat Neutral Transfert Coefficients 1310 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)1311 1562 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 1563 1312 1564 ! Atmospheric and Surface Variables 1313 1565 zst(:,:) = sst_m(:,:) + rt0 ! convert SST from Celcius to Kelvin 1314 zqo_sat(:,:) = 0.98_wp * q_sat( zst(:,:) , sf(jp_slp)%fnow(:,:,1) ) ! saturation humidity over ocean [kg/kg] 1315 zqi_sat(:,:) = 0.98_wp * q_sat( ztm_su(:,:), sf(jp_slp)%fnow(:,:,1) ) ! saturation humidity over ice [kg/kg] 1316 ! 1317 DO jj = 2, jpjm1 ! reduced loop is necessary for reproducibility 1318 DO ji = fs_2, fs_jpim1 1319 ! Virtual potential temperature [K] 1320 zthetav_os = zst(ji,jj) * ( 1._wp + rctv0 * zqo_sat(ji,jj) ) ! over ocean 1321 zthetav_is = ztm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) ) ! ocean ice 1322 zthetav_zu = t_zu (ji,jj) * ( 1._wp + rctv0 * q_zu(ji,jj) ) ! at zu 1323 1324 ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead) 1325 zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj) )**2 ! over ocean 1326 zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2 ! over ice 1327 1328 ! Momentum and Heat Neutral Transfert Coefficients 1329 zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta ! Eq. 40 1330 zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) ) ! Eq. 53 1331 1332 ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead) 1333 z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water 1334 z0i = z0_skin_ice ! over ice (cf Lupkes email for details) 1335 IF( zrib_o <= 0._wp ) THEN 1336 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 1337 zfhw = ( 1._wp + ( zbetah * ( zthetav_os - zthetav_zu )**r1_3 / ( Chn_oce(ji,jj) * MAX(0.01, wndm(ji,jj)) ) & ! Eq. 26 1338 & )**zgamma )**z1_gamma 1339 ELSE 1340 zfmw = 1._wp / ( 1._wp + zam * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 12 1341 zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 28 1342 ENDIF 1343 1344 IF( zrib_i <= 0._wp ) THEN 1345 zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 9 1346 zfhi = 1._wp - zah * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 25 1347 ELSE 1348 zfmi = 1._wp / ( 1._wp + zam * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 11 1349 zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 27 1350 ENDIF 1351 1352 ! Momentum Transfert Coefficients (Eq. 38) 1353 Cd(ji,jj) = zCdn_skin_ice * zfmi + & 1354 & 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) ) 1355 1356 ! Heat Transfert Coefficients (Eq. 49) 1357 Ch(ji,jj) = zChn_skin_ice * zfhi + & 1358 & 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) ) 1359 ! 1360 END DO 1361 END DO 1362 CALL lbc_lnk_multi( 'sbcblk', Cd, 'T', 1., Ch, 'T', 1. ) 1566 zqo_sat(:,:) = rdct_qsat_salt * q_sat( zst(:,:) , pslp(:,:) ) ! saturation humidity over ocean [kg/kg] 1567 zqi_sat(:,:) = q_sat( ptm_su(:,:), pslp(:,:) ) ! saturation humidity over ice [kg/kg] 1568 ! 1569 DO_2D( 0, 0, 0, 0 ) 1570 ! Virtual potential temperature [K] 1571 zthetav_os = zst(ji,jj) * ( 1._wp + rctv0 * zqo_sat(ji,jj) ) ! over ocean 1572 zthetav_is = ptm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) ) ! ocean ice 1573 zthetav_zu = t_zu (ji,jj) * ( 1._wp + rctv0 * q_zu(ji,jj) ) ! at zu 1574 1575 ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead) 1576 zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj) )**2 ! over ocean 1577 zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2 ! over ice 1578 1579 ! Momentum and Heat Neutral Transfert Coefficients 1580 zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta ! Eq. 40 1581 zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) ) ! Eq. 53 1582 1583 ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead ?) 1584 z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water 1585 z0i = z0_skin_ice ! over ice 1586 IF( zrib_o <= 0._wp ) THEN 1587 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 1588 zfhw = ( 1._wp + ( zbetah * ( zthetav_os - zthetav_zu )**r1_3 / ( Chn_oce(ji,jj) * MAX(0.01, wndm(ji,jj)) ) & ! Eq. 26 1589 & )**zgamma )**z1_gamma 1590 ELSE 1591 zfmw = 1._wp / ( 1._wp + zam * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 12 1592 zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 28 1593 ENDIF 1594 1595 IF( zrib_i <= 0._wp ) THEN 1596 zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 9 1597 zfhi = 1._wp - zah * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 25 1598 ELSE 1599 zfmi = 1._wp / ( 1._wp + zam * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 11 1600 zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 27 1601 ENDIF 1602 1603 ! Momentum Transfert Coefficients (Eq. 38) 1604 pcd(ji,jj) = zCdn_skin_ice * zfmi + & 1605 & 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) ) 1606 1607 ! Heat Transfert Coefficients (Eq. 49) 1608 pch(ji,jj) = zChn_skin_ice * zfhi + & 1609 & 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) ) 1610 ! 1611 END_2D 1612 CALL lbc_lnk_multi( 'sbcblk', pcd, 'T', 1.0_wp, pch, 'T', 1.0_wp ) 1363 1613 ! 1364 1614 END SUBROUTINE Cdn10_Lupkes2015
Note: See TracChangeset
for help on using the changeset viewer.
