- Timestamp:
- 2019-12-10T15:44:23+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/SBC/sbcblk.F90
r12109 r12154 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 :: jpfld ! maximum number of files to read 77 INTEGER , PUBLIC, PARAMETER :: jp_wndi = 1 ! index of 10m wind velocity (i-component) (m/s) at T-point 78 INTEGER , PUBLIC, PARAMETER :: jp_wndj = 2 ! index of 10m wind velocity (j-component) (m/s) at T-point 79 INTEGER , PUBLIC, PARAMETER :: jp_tair = 3 ! index of 10m air temperature (Kelvin) 80 INTEGER , PUBLIC, PARAMETER :: jp_humi = 4 ! index of specific humidity ( % ) 81 INTEGER , PUBLIC, PARAMETER :: jp_qsr = 5 ! index of solar heat (W/m2) 82 INTEGER , PUBLIC, PARAMETER :: jp_qlw = 6 ! index of Long wave (W/m2) 83 INTEGER , PUBLIC, PARAMETER :: jp_prec = 7 ! index of total precipitation (rain+snow) (Kg/m2/s) 84 INTEGER , PUBLIC, PARAMETER :: jp_snow = 8 ! index of snow (solid prcipitation) (kg/m2/s) 85 INTEGER , PUBLIC, PARAMETER :: jp_slp = 9 ! index of sea level pressure (Pa) 86 INTEGER , PUBLIC, PARAMETER :: jp_hpgi =10 ! index of ABL geostrophic wind or hpg (i-component) (m/s) at T-point 87 INTEGER , PUBLIC, PARAMETER :: jp_hpgj =11 ! index of ABL geostrophic wind or hpg (j-component) (m/s) at T-point 88 89 TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) :: sf ! structure of input atmospheric fields (file informations, fields read) 90 103 91 ! !!* Namelist namsbc_blk : bulk parameters 104 92 LOGICAL :: ln_NCAR ! "NCAR" algorithm (Large and Yeager 2008) 105 93 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)94 LOGICAL :: ln_COARE_3p6 ! "COARE 3.6" algorithm (Edson et al. 2013) 95 LOGICAL :: ln_ECMWF ! "ECMWF" algorithm (IFS cycle 45r1) 108 96 ! 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) 97 LOGICAL :: ln_Cd_L12 ! ice-atm drag = F( ice concentration ) (Lupkes et al. JGR2012) 98 LOGICAL :: ln_Cd_L15 ! ice-atm drag = F( ice concentration, atmospheric stability ) (Lupkes et al. JGR2015) 118 99 ! 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 100 REAL(wp) :: rn_pfac ! multiplication factor for precipitation 101 REAL(wp), PUBLIC :: rn_efac ! multiplication factor for evaporation 102 REAL(wp), PUBLIC :: rn_vfac ! multiplication factor for ice/ocean velocity in the calculation of wind stress 103 REAL(wp) :: rn_zqt ! z(q,t) : height of humidity and temperature measurements 104 REAL(wp) :: rn_zu ! z(u) : height of wind measurements 105 ! 106 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: Cd_ice , Ch_ice , Ce_ice ! transfert coefficients over ice 107 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: Cdn_oce, Chn_oce, Cen_oce ! neutral coeffs over ocean (L15 bulk scheme) 108 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: t_zu, q_zu ! air temp. and spec. hum. at wind speed height (L15 bulk scheme) 109 110 LOGICAL :: ln_skin_cs ! use the cool-skin (only available in ECMWF and COARE algorithms) !LB 111 LOGICAL :: ln_skin_wl ! use the warm-layer parameterization (only available in ECMWF and COARE algorithms) !LB 112 LOGICAL :: ln_humi_sph ! humidity read in files ("sn_humi") is specific humidity [kg/kg] if .true. !LB 113 LOGICAL :: ln_humi_dpt ! humidity read in files ("sn_humi") is dew-point temperature [K] if .true. !LB 114 LOGICAL :: ln_humi_rlh ! humidity read in files ("sn_humi") is relative humidity [%] if .true. !LB 115 ! 116 INTEGER :: nhumi ! choice of the bulk algorithm 117 ! ! associated indices: 118 INTEGER, PARAMETER :: np_humi_sph = 1 119 INTEGER, PARAMETER :: np_humi_dpt = 2 120 INTEGER, PARAMETER :: np_humi_rlh = 3 125 121 126 122 INTEGER :: nblk ! choice of the bulk algorithm … … 128 124 INTEGER, PARAMETER :: np_NCAR = 1 ! "NCAR" algorithm (Large and Yeager 2008) 129 125 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)126 INTEGER, PARAMETER :: np_COARE_3p6 = 3 ! "COARE 3.6" algorithm (Edson et al. 2013) 127 INTEGER, PARAMETER :: np_ECMWF = 4 ! "ECMWF" algorithm (IFS cycle 45r1) 132 128 133 129 !! * Substitutions … … 144 140 !! *** ROUTINE sbc_blk_alloc *** 145 141 !!------------------------------------------------------------------- 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 ) 142 ALLOCATE( t_zu(jpi,jpj) , q_zu(jpi,jpj) , & 143 & Cdn_oce(jpi,jpj), Chn_oce(jpi,jpj), Cen_oce(jpi,jpj), & 144 & Cd_ice (jpi,jpj), Ch_ice (jpi,jpj), Ce_ice (jpi,jpj), STAT=sbc_blk_alloc ) 148 145 ! 149 146 CALL mpp_sum ( 'sbcblk', sbc_blk_alloc ) … … 158 155 !! ** Purpose : choose and initialize a bulk formulae formulation 159 156 !! 160 !! ** Method : 157 !! ** Method : 161 158 !! 162 159 !!---------------------------------------------------------------------- 163 INTEGER :: ifpr, jfld! dummy loop indice and argument160 INTEGER :: jfpr ! dummy loop indice and argument 164 161 INTEGER :: ios, ierror, ioptio ! Local integer 165 162 !! 166 163 CHARACTER(len=100) :: cn_dir ! Root directory for location of atmospheric forcing files 167 TYPE(FLD_N), DIMENSION(jpfld) :: slf_i! array of namelist informations on the fields to read164 TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) :: slf_i ! array of namelist informations on the fields to read 168 165 TYPE(FLD_N) :: sn_wndi, sn_wndj, sn_humi, sn_qsr ! informations about the fields to be read 169 166 TYPE(FLD_N) :: sn_qlw , sn_tair, sn_prec, sn_snow ! " " 170 TYPE(FLD_N) :: sn_slp , sn_ tdif! " "167 TYPE(FLD_N) :: sn_slp , sn_hpgi, sn_hpgj ! " " 171 168 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 169 & sn_tair, sn_prec, sn_snow, sn_slp, sn_hpgi, sn_hpgj, & 170 & ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF, & ! bulk algorithm 171 & cn_dir , rn_zqt, rn_zu, & 172 & rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15, & 173 & ln_skin_cs, ln_skin_wl, ln_humi_sph, ln_humi_dpt, ln_humi_rlh ! cool-skin / warm-layer !LB 176 174 !!--------------------------------------------------------------------- 177 175 ! … … 179 177 IF( sbc_blk_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' ) 180 178 ! 181 ! !** read bulk namelist 179 ! !** read bulk namelist 182 180 REWIND( numnam_ref ) !* Namelist namsbc_blk in reference namelist : bulk parameters 183 181 READ ( numnam_ref, namsbc_blk, IOSTAT = ios, ERR = 901) … … 192 190 ! !** initialization of the chosen bulk formulae (+ check) 193 191 ! !* 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 ! 192 ioptio = 0 193 IF( ln_NCAR ) THEN 194 nblk = np_NCAR ; ioptio = ioptio + 1 195 ENDIF 196 IF( ln_COARE_3p0 ) THEN 197 nblk = np_COARE_3p0 ; ioptio = ioptio + 1 198 ENDIF 199 IF( ln_COARE_3p6 ) THEN 200 nblk = np_COARE_3p6 ; ioptio = ioptio + 1 201 ENDIF 202 IF( ln_ECMWF ) THEN 203 nblk = np_ECMWF ; ioptio = ioptio + 1 204 ENDIF 200 205 IF( ioptio /= 1 ) CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' ) 206 207 ! !** initialization of the cool-skin / warm-layer parametrization 208 IF( ln_skin_cs .OR. ln_skin_wl ) THEN 209 !! Some namelist sanity tests: 210 IF( ln_NCAR ) & 211 & CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with NCAR algorithm' ) 212 IF( nn_fsbc /= 1 ) & 213 & CALL ctl_stop( 'sbc_blk_init: Please set "nn_fsbc" to 1 when using cool-skin/warm-layer param.') 214 END IF 215 216 IF( ln_skin_wl ) THEN 217 !! Check if the frequency of downwelling solar flux input makes sense and if ln_dm2dc=T if it is daily! 218 IF( (sn_qsr%freqh < 0.).OR.(sn_qsr%freqh > 24.) ) & 219 & CALL ctl_stop( 'sbc_blk_init: Warm-layer param. (ln_skin_wl) not compatible with freq. of solar flux > daily' ) 220 IF( (sn_qsr%freqh == 24.).AND.(.NOT. ln_dm2dc) ) & 221 & CALL ctl_stop( 'sbc_blk_init: Please set ln_dm2dc=T for warm-layer param. (ln_skin_wl) to work properly' ) 222 END IF 223 224 ioptio = 0 225 IF( ln_humi_sph ) THEN 226 nhumi = np_humi_sph ; ioptio = ioptio + 1 227 ENDIF 228 IF( ln_humi_dpt ) THEN 229 nhumi = np_humi_dpt ; ioptio = ioptio + 1 230 ENDIF 231 IF( ln_humi_rlh ) THEN 232 nhumi = np_humi_rlh ; ioptio = ioptio + 1 233 ENDIF 234 IF( ioptio /= 1 ) CALL ctl_stop( 'sbc_blk_init: Choose one and only one type of air humidity' ) 201 235 ! 202 236 IF( ln_dm2dc ) THEN !* check: diurnal cycle on Qsr 203 237 IF( sn_qsr%freqh /= 24. ) CALL ctl_stop( 'sbc_blk_init: ln_dm2dc=T only with daily short-wave input' ) 204 IF( sn_qsr%ln_tint ) THEN 238 IF( sn_qsr%ln_tint ) THEN 205 239 CALL ctl_warn( 'sbc_blk_init: ln_dm2dc=T daily qsr time interpolation done by sbcdcy module', & 206 240 & ' ==> We force time interpolation = .false. for qsr' ) … … 210 244 ! !* set the bulk structure 211 245 ! !- store namelist information in an array 246 IF( ln_blk ) jpfld = 9 247 IF( ln_abl ) jpfld = 11 248 ALLOCATE( slf_i(jpfld) ) 249 ! 212 250 slf_i(jp_wndi) = sn_wndi ; slf_i(jp_wndj) = sn_wndj 213 251 slf_i(jp_qsr ) = sn_qsr ; slf_i(jp_qlw ) = sn_qlw 214 252 slf_i(jp_tair) = sn_tair ; slf_i(jp_humi) = sn_humi 215 253 slf_i(jp_prec) = sn_prec ; slf_i(jp_snow) = sn_snow 216 slf_i(jp_slp ) = sn_slp ; slf_i(jp_tdif) = sn_tdif217 !218 lhftau = ln_taudif !- add an extra field if HF stress is used219 jfld = jpfld - COUNT( (/.NOT.lhftau/) )254 slf_i(jp_slp ) = sn_slp 255 IF( ln_abl ) THEN 256 slf_i(jp_hpgi) = sn_hpgi ; slf_i(jp_hpgj) = sn_hpgj 257 END IF 220 258 ! 221 259 ! !- allocate the bulk structure 222 ALLOCATE( sf(j fld), STAT=ierror )260 ALLOCATE( sf(jpfld), STAT=ierror ) 223 261 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)%freqh > 0. .AND. MOD( NINT(3600. * slf_i(ifpr)%freqh), nn_fsbc * NINT(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 262 ! 263 DO jfpr= 1, jpfld 264 ! 265 IF( TRIM(sf(jfpr)%clrootname) == 'NOT USED' ) THEN !-- not used field --! (only now allocated and set to zero) 266 ALLOCATE( sf(jfpr)%fnow(jpi,jpj,1) ) 267 sf(jfpr)%fnow(:,:,1) = 0._wp 268 ELSE !-- used field --! 269 IF( ln_abl .AND. & 270 & ( jfpr == jp_wndi .OR. jfpr == jp_wndj .OR. jfpr == jp_humi .OR. & 271 & jfpr == jp_hpgi .OR. jfpr == jp_hpgj .OR. jfpr == jp_tair ) ) THEN ! ABL: some fields are 3D input 272 ALLOCATE( sf(jfpr)%fnow(jpi,jpj,jpka) ) 273 IF( slf_i(jfpr)%ln_tint ) ALLOCATE( sf(jfpr)%fdta(jpi,jpj,jpka,2) ) 274 ELSE ! others or Bulk fields are 2D fiels 275 ALLOCATE( sf(jfpr)%fnow(jpi,jpj,1) ) 276 IF( slf_i(jfpr)%ln_tint ) ALLOCATE( sf(jfpr)%fdta(jpi,jpj,1,2) ) 277 ENDIF 278 ! 279 IF( slf_i(jfpr)%freqh > 0. .AND. MOD( NINT(3600. * slf_i(jfpr)%freqh), nn_fsbc * NINT(rdt) ) /= 0 ) & 280 & CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rdt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.', & 281 & ' This is not ideal. You should consider changing either rdt or nn_fsbc value...' ) 282 ENDIF 231 283 END DO 232 284 ! !- fill the bulk structure with namelist informations 233 285 CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' ) 234 286 ! 235 IF 236 !Activated wave module but neither drag nor stokes drift activated237 IF 287 IF( ln_wave ) THEN 288 !Activated wave module but neither drag nor stokes drift activated 289 IF( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) ) THEN 238 290 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 291 !drag coefficient read from wave model definable only with mfs bulk formulae and core 292 ELSEIF(ln_cdgw .AND. .NOT. ln_NCAR ) THEN 293 CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR and CORE bulk formulae') 294 ELSEIF(ln_stcor .AND. .NOT. ln_sdw) THEN 295 CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T') 244 296 ENDIF 245 297 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 ! 298 IF( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) & 299 & CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ', & 300 & 'with drag coefficient (ln_cdgw =T) ' , & 301 & 'or Stokes Drift (ln_sdw=T) ' , & 302 & 'or ocean stress modification due to waves (ln_tauwoc=T) ', & 303 & 'or Stokes-Coriolis term (ln_stcori=T)' ) 304 ENDIF 305 ! 306 IF( ln_abl ) THEN ! ABL: read 3D fields for wind, temperature, humidity and pressure gradient 307 rn_zqt = ght_abl(2) ! set the bulk altitude to ABL first level 308 rn_zu = ght_abl(2) 309 IF(lwp) WRITE(numout,*) 310 IF(lwp) WRITE(numout,*) ' ABL formulation: overwrite rn_zqt & rn_zu with ABL first level altitude' 311 ENDIF 312 ! 313 ! set transfer coefficients to default sea-ice values 314 Cd_ice(:,:) = rCd_ice 315 Ch_ice(:,:) = rCd_ice 316 Ce_ice(:,:) = rCd_ice 317 ! 255 318 IF(lwp) THEN !** Control print 256 319 ! 257 WRITE(numout,*) !* namelist 320 WRITE(numout,*) !* namelist 258 321 WRITE(numout,*) ' Namelist namsbc_blk (other than data information):' 259 322 WRITE(numout,*) ' "NCAR" algorithm (Large and Yeager 2008) ln_NCAR = ', ln_NCAR 260 323 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 324 WRITE(numout,*) ' "COARE 3.6" algorithm (Fairall 2018 + Edson al 2013)ln_COARE_3p6 = ', ln_COARE_3p6 325 WRITE(numout,*) ' "ECMWF" algorithm (IFS cycle 45r1) ln_ECMWF = ', ln_ECMWF 264 326 WRITE(numout,*) ' Air temperature and humidity reference height (m) rn_zqt = ', rn_zqt 265 327 WRITE(numout,*) ' Wind vector reference height (m) rn_zu = ', rn_zu … … 275 337 CASE( np_NCAR ) ; WRITE(numout,*) ' ==>>> "NCAR" algorithm (Large and Yeager 2008)' 276 338 CASE( np_COARE_3p0 ) ; WRITE(numout,*) ' ==>>> "COARE 3.0" algorithm (Fairall et al. 2003)' 277 CASE( np_COARE_3p 5 ) ; WRITE(numout,*) ' ==>>> "COARE 3.5" algorithm (Edson et al. 2013)'278 CASE( np_ECMWF ) ; WRITE(numout,*) ' ==>>> "ECMWF" algorithm (IFS cycle 31)'339 CASE( np_COARE_3p6 ) ; WRITE(numout,*) ' ==>>> "COARE 3.6" algorithm (Fairall 2018+Edson et al. 2013)' 340 CASE( np_ECMWF ) ; WRITE(numout,*) ' ==>>> "ECMWF" algorithm (IFS cycle 45r1)' 279 341 END SELECT 280 342 ! 343 WRITE(numout,*) 344 WRITE(numout,*) ' use cool-skin parameterization (SSST) ln_skin_cs = ', ln_skin_cs 345 WRITE(numout,*) ' use warm-layer parameterization (SSST) ln_skin_wl = ', ln_skin_wl 346 ! 347 WRITE(numout,*) 348 SELECT CASE( nhumi ) !* Print the choice of air humidity 349 CASE( np_humi_sph ) ; WRITE(numout,*) ' ==>>> air humidity is SPECIFIC HUMIDITY [kg/kg]' 350 CASE( np_humi_dpt ) ; WRITE(numout,*) ' ==>>> air humidity is DEW-POINT TEMPERATURE [K]' 351 CASE( np_humi_rlh ) ; WRITE(numout,*) ' ==>>> air humidity is RELATIVE HUMIDITY [%]' 352 END SELECT 353 ! 281 354 ENDIF 282 355 ! … … 291 364 !! (momentum, heat, freshwater and runoff) 292 365 !! 293 !! ** Method : (1) READ each fluxes in NetCDF files: 294 !! the 10m wind velocity (i-component) (m/s) at T-point 295 !! the 10m wind velocity (j-component) (m/s) at T-point 296 !! the 10m or 2m specific humidity ( % ) 297 !! the solar heat (W/m2) 298 !! the Long wave (W/m2) 299 !! the 10m or 2m air temperature (Kelvin) 300 !! the total precipitation (rain+snow) (Kg/m2/s) 301 !! the snow (solid prcipitation) (kg/m2/s) 302 !! the tau diff associated to HF tau (N/m2) at T-point (ln_taudif=T) 303 !! (2) CALL blk_oce 366 !! ** Method : 367 !! (1) READ each fluxes in NetCDF files: 368 !! the wind velocity (i-component) at z=rn_zu (m/s) at T-point 369 !! the wind velocity (j-component) at z=rn_zu (m/s) at T-point 370 !! the specific humidity at z=rn_zqt (kg/kg) 371 !! the air temperature at z=rn_zqt (Kelvin) 372 !! the solar heat (W/m2) 373 !! the Long wave (W/m2) 374 !! the total precipitation (rain+snow) (Kg/m2/s) 375 !! the snow (solid precipitation) (kg/m2/s) 376 !! ABL dynamical forcing (i/j-components of either hpg or geostrophic winds) 377 !! (2) CALL blk_oce_1 and blk_oce_2 304 378 !! 305 379 !! C A U T I O N : never mask the surface stress fields … … 318 392 !!---------------------------------------------------------------------- 319 393 INTEGER, INTENT(in) :: kt ! ocean time step 320 !!--------------------------------------------------------------------- 394 !!---------------------------------------------------------------------- 395 REAL(wp), DIMENSION(jpi,jpj) :: zssq, zcd_du, zsen, zevp 396 REAL(wp) :: ztmp 397 !!---------------------------------------------------------------------- 321 398 ! 322 399 CALL fld_read( kt, nn_fsbc, sf ) ! input fields provided at the current time-step 323 ! 400 401 ! Sanity/consistence test on humidity at first time step to detect potential screw-up: 402 IF( kt == nit000 ) THEN 403 WRITE(numout,*) '' 404 #if defined key_agrif 405 WRITE(numout,*) ' === AGRIF => Sanity/consistence test on air humidity SKIPPED! :( ===' 406 #else 407 ztmp = SUM(tmask(:,:,1)) ! number of ocean points on local proc domain 408 IF( ztmp > 8._wp ) THEN ! test only on proc domains with at least 8 ocean points! 409 ztmp = SUM(sf(jp_humi)%fnow(:,:,1)*tmask(:,:,1))/ztmp ! mean humidity over ocean on proc 410 SELECT CASE( nhumi ) 411 CASE( np_humi_sph ) ! specific humidity => expect: 0. <= something < 0.065 [kg/kg] (0.061 is saturation at 45degC !!!) 412 IF( (ztmp < 0._wp) .OR. (ztmp > 0.065) ) ztmp = -1._wp 413 CASE( np_humi_dpt ) ! dew-point temperature => expect: 110. <= something < 320. [K] 414 IF( (ztmp < 110._wp).OR.(ztmp > 320._wp) ) ztmp = -1._wp 415 CASE( np_humi_rlh ) ! relative humidity => expect: 0. <= something < 100. [%] 416 IF( (ztmp < 0._wp) .OR.(ztmp > 100._wp) ) ztmp = -1._wp 417 END SELECT 418 IF(ztmp < 0._wp) THEN 419 WRITE(numout,'(" Mean humidity value found on proc #",i5.5," is: ",f)'), narea, ztmp 420 CALL ctl_stop( 'STOP', 'Something is wrong with air humidity!!!', & 421 & ' ==> check the unit in your input files' , & 422 & ' ==> check consistence of namelist choice: specific? relative? dew-point?', & 423 & ' ==> ln_humi_sph -> [kg/kg] | ln_humi_rlh -> [%] | ln_humi_dpt -> [K] !!!' ) 424 END IF 425 END IF 426 WRITE(numout,*) ' === Sanity/consistence test on air humidity sucessfuly passed! ===' 427 #endif 428 WRITE(numout,*) '' 429 END IF !IF( kt == nit000 ) 324 430 ! ! compute the surface ocean fluxes using bulk formulea 325 IF( MOD( kt - 1, nn_fsbc ) == 0 ) CALL blk_oce( kt, sf, sst_m, ssu_m, ssv_m ) 326 431 IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 432 CALL blk_oce_1( kt, sf(jp_wndi)%fnow(:,:,1), sf(jp_wndj)%fnow(:,:,1), & ! <<= in 433 & sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), & ! <<= in 434 & sf(jp_slp )%fnow(:,:,1), sst_m, ssu_m, ssv_m, & ! <<= in 435 & sf(jp_qsr )%fnow(:,:,1), sf(jp_qlw )%fnow(:,:,1), & ! <<= in (wl/cs) 436 & zssq, zcd_du, zsen, zevp ) ! =>> out 437 438 CALL blk_oce_2( sf(jp_tair)%fnow(:,:,1), sf(jp_qsr )%fnow(:,:,1), & ! <<= in 439 & sf(jp_qlw )%fnow(:,:,1), sf(jp_prec)%fnow(:,:,1), & ! <<= in 440 & sf(jp_snow)%fnow(:,:,1), sst_m, & ! <<= in 441 & zsen, zevp ) ! <=> in out 442 ENDIF 443 ! 327 444 #if defined key_cice 328 445 IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 329 446 qlw_ice(:,:,1) = sf(jp_qlw )%fnow(:,:,1) 330 IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) 331 ELSE ; qsr_ice(:,:,1) = sf(jp_qsr)%fnow(:,:,1) 332 ENDIF 447 IF( ln_dm2dc ) THEN 448 qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) 449 ELSE 450 qsr_ice(:,:,1) = sf(jp_qsr)%fnow(:,:,1) 451 ENDIF 333 452 tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) 334 qatm_ice(:,:) = sf(jp_humi)%fnow(:,:,1) 453 454 SELECT CASE( nhumi ) 455 CASE( np_humi_sph ) 456 qatm_ice(:,:) = sf(jp_humi)%fnow(:,:,1) 457 CASE( np_humi_dpt ) 458 qatm_ice(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 459 CASE( np_humi_rlh ) 460 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 461 END SELECT 462 335 463 tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac 336 464 sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac … … 343 471 344 472 345 SUBROUTINE blk_oce( kt, sf, pst, pu, pv ) 346 !!--------------------------------------------------------------------- 347 !! *** ROUTINE blk_oce *** 348 !! 349 !! ** Purpose : provide the momentum, heat and freshwater fluxes at 350 !! the ocean surface at each time step 351 !! 352 !! ** Method : bulk formulea for the ocean using atmospheric 353 !! fields read in sbc_read 473 SUBROUTINE blk_oce_1( kt, pwndi, pwndj , ptair, phumi, & ! inp 474 & pslp , pst , pu , pv, & ! inp 475 & pqsr , pqlw , & ! inp 476 & pssq , pcd_du, psen , pevp ) ! out 477 !!--------------------------------------------------------------------- 478 !! *** ROUTINE blk_oce_1 *** 479 !! 480 !! ** Purpose : if ln_blk=T, computes surface momentum, heat and freshwater fluxes 481 !! if ln_abl=T, computes Cd x |U|, Ch x |U|, Ce x |U| for ABL integration 482 !! 483 !! ** Method : bulk formulae using atmospheric fields from : 484 !! if ln_blk=T, atmospheric fields read in sbc_read 485 !! if ln_abl=T, the ABL model at previous time-step 486 !! 487 !! ** Outputs : - pssq : surface humidity used to compute latent heat flux (kg/kg) 488 !! - pcd_du : Cd x |dU| at T-points (m/s) 489 !! - psen : Ch x |dU| at T-points (m/s) 490 !! - pevp : Ce x |dU| at T-points (m/s) 491 !!--------------------------------------------------------------------- 492 INTEGER , INTENT(in ) :: kt ! time step index 493 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pwndi ! atmospheric wind at U-point [m/s] 494 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pwndj ! atmospheric wind at V-point [m/s] 495 REAL(wp), INTENT(in ), DIMENSION(:,:) :: phumi ! specific humidity at T-points [kg/kg] 496 REAL(wp), INTENT(in ), DIMENSION(:,:) :: ptair ! potential temperature at T-points [Kelvin] 497 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pslp ! sea-level pressure [Pa] 498 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pst ! surface temperature [Celcius] 499 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pu ! surface current at U-point (i-component) [m/s] 500 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pv ! surface current at V-point (j-component) [m/s] 501 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pqsr ! 502 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pqlw ! 503 REAL(wp), INTENT( out), DIMENSION(:,:) :: pssq ! specific humidity at pst [kg/kg] 504 REAL(wp), INTENT( out), DIMENSION(:,:) :: pcd_du ! Cd x |dU| at T-points [m/s] 505 REAL(wp), INTENT( out), DIMENSION(:,:) :: psen ! Ch x |dU| at T-points [m/s] 506 REAL(wp), INTENT( out), DIMENSION(:,:) :: pevp ! Ce x |dU| at T-points [m/s] 507 ! 508 INTEGER :: ji, jj ! dummy loop indices 509 REAL(wp) :: zztmp ! local variable 510 REAL(wp), DIMENSION(jpi,jpj) :: zwnd_i, zwnd_j ! wind speed components at T-point 511 REAL(wp), DIMENSION(jpi,jpj) :: zst ! surface temperature in Kelvin 512 REAL(wp), DIMENSION(jpi,jpj) :: zU_zu ! bulk wind speed at height zu [m/s] 513 REAL(wp), DIMENSION(jpi,jpj) :: ztpot ! potential temperature of air at z=rn_zqt [K] 514 REAL(wp), DIMENSION(jpi,jpj) :: zqair ! specific humidity of air at z=rn_zqt [kg/kg] 515 REAL(wp), DIMENSION(jpi,jpj) :: zcd_oce ! momentum transfert coefficient over ocean 516 REAL(wp), DIMENSION(jpi,jpj) :: zch_oce ! sensible heat transfert coefficient over ocean 517 REAL(wp), DIMENSION(jpi,jpj) :: zce_oce ! latent heat transfert coefficient over ocean 518 REAL(wp), DIMENSION(jpi,jpj) :: zqla ! latent heat flux 519 REAL(wp), DIMENSION(jpi,jpj) :: zztmp1, zztmp2 520 !!--------------------------------------------------------------------- 521 ! 522 ! local scalars ( place there for vector optimisation purposes) 523 zst(:,:) = pst(:,:) + rt0 ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K) 524 525 ! ----------------------------------------------------------------------------- ! 526 ! 0 Wind components and module at T-point relative to the moving ocean ! 527 ! ----------------------------------------------------------------------------- ! 528 529 ! ... components ( U10m - U_oce ) at T-point (unmasked) 530 #if defined key_cyclone 531 zwnd_i(:,:) = 0._wp 532 zwnd_j(:,:) = 0._wp 533 CALL wnd_cyc( kt, zwnd_i, zwnd_j ) ! add analytical tropical cyclone (Vincent et al. JGR 2012) 534 DO jj = 2, jpjm1 535 DO ji = fs_2, fs_jpim1 ! vect. opt. 536 pwndi(ji,jj) = pwndi(ji,jj) + zwnd_i(ji,jj) 537 pwndj(ji,jj) = pwndj(ji,jj) + zwnd_j(ji,jj) 538 END DO 539 END DO 540 #endif 541 DO jj = 2, jpjm1 542 DO ji = fs_2, fs_jpim1 ! vect. opt. 543 zwnd_i(ji,jj) = ( pwndi(ji,jj) - rn_vfac * 0.5 * ( pu(ji-1,jj ) + pu(ji,jj) ) ) 544 zwnd_j(ji,jj) = ( pwndj(ji,jj) - rn_vfac * 0.5 * ( pv(ji ,jj-1) + pv(ji,jj) ) ) 545 END DO 546 END DO 547 CALL lbc_lnk_multi( 'sbcblk', zwnd_i, 'T', -1., zwnd_j, 'T', -1. ) 548 ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked) 549 wndm(:,:) = SQRT( zwnd_i(:,:) * zwnd_i(:,:) & 550 & + zwnd_j(:,:) * zwnd_j(:,:) ) * tmask(:,:,1) 551 552 ! ----------------------------------------------------------------------------- ! 553 ! I Solar FLUX ! 554 ! ----------------------------------------------------------------------------- ! 555 556 ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle ! Short Wave 557 zztmp = 1. - albo 558 IF( ln_dm2dc ) THEN 559 qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1) 560 ELSE 561 qsr(:,:) = zztmp * sf(jp_qsr)%fnow(:,:,1) * tmask(:,:,1) 562 ENDIF 563 564 565 ! ----------------------------------------------------------------------------- ! 566 ! II Turbulent FLUXES ! 567 ! ----------------------------------------------------------------------------- ! 568 569 ! specific humidity at SST 570 pssq(:,:) = rdct_qsat_salt * q_sat( zst(:,:), pslp(:,:) ) 571 572 IF( ln_skin_cs .OR. ln_skin_wl ) THEN 573 zztmp1(:,:) = zst(:,:) 574 zztmp2(:,:) = pssq(:,:) 575 ENDIF 576 577 ! specific humidity of air at "rn_zqt" m above the sea 578 SELECT CASE( nhumi ) 579 CASE( np_humi_sph ) 580 zqair(:,:) = phumi(:,:) ! what we read in file is already a spec. humidity! 581 CASE( np_humi_dpt ) 582 !IF(lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of d_air and slp !' !LBrm 583 zqair(:,:) = q_sat( phumi(:,:), pslp(:,:) ) 584 CASE( np_humi_rlh ) 585 !IF(lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of RH, t_air and slp !' !LBrm 586 zqair(:,:) = q_air_rh( 0.01_wp*phumi(:,:), ptair(:,:), pslp(:,:) ) !LB: 0.01 => RH is % percent in file 587 END SELECT 588 ! 589 ! potential temperature of air at "rn_zqt" m above the sea 590 IF( ln_abl ) THEN 591 ztpot = ptair(:,:) 592 ELSE 593 ! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate 594 ! (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2 595 ! (since reanalysis products provide T at z, not theta !) 596 !#LB: because AGRIF hates functions that return something else than a scalar, need to 597 ! use scalar version of gamma_moist() ... 598 DO jj = 1, jpj 599 DO ji = 1, jpi 600 ztpot(ji,jj) = ptair(ji,jj) + gamma_moist( ptair(ji,jj), zqair(ji,jj) ) * rn_zqt 601 END DO 602 END DO 603 ENDIF 604 605 606 607 !! Time to call the user-selected bulk parameterization for 608 !! == transfer coefficients ==! Cd, Ch, Ce at T-point, and more... 609 SELECT CASE( nblk ) 610 611 CASE( np_NCAR ) 612 CALL turb_ncar ( rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, & 613 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 614 615 CASE( np_COARE_3p0 ) 616 CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 617 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 618 & Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 619 620 CASE( np_COARE_3p6 ) 621 CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 622 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 623 & Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 624 625 CASE( np_ECMWF ) 626 CALL turb_ecmwf ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 627 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 628 & Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 629 630 CASE DEFAULT 631 CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 632 633 END SELECT 634 635 IF( ln_skin_cs .OR. ln_skin_wl ) THEN 636 !! In the presence of sea-ice we forget about the cool-skin/warm-layer update of zst and pssq: 637 WHERE ( fr_i < 0.001_wp ) 638 ! zst and pssq have been updated by cool-skin/warm-layer scheme and we keep it!!! 639 zst(:,:) = zst(:,:)*tmask(:,:,1) 640 pssq(:,:) = pssq(:,:)*tmask(:,:,1) 641 ELSEWHERE 642 ! we forget about the update... 643 zst(:,:) = zztmp1(:,:) !#LB: using what we backed up before skin-algo 644 pssq(:,:) = zztmp2(:,:) !#LB: " " " 645 END WHERE 646 END IF 647 648 !! CALL iom_put( "Cd_oce", zcd_oce) ! output value of pure ocean-atm. transfer coef. 649 !! CALL iom_put( "Ch_oce", zch_oce) ! output value of pure ocean-atm. transfer coef. 650 651 IF( ABS(rn_zu - rn_zqt) < 0.1_wp ) THEN 652 !! If zu == zt, then ensuring once for all that: 653 t_zu(:,:) = ztpot(:,:) 654 q_zu(:,:) = zqair(:,:) 655 ENDIF 656 657 658 ! Turbulent fluxes over ocean => BULK_FORMULA @ sbcblk_phy.F90 659 ! ------------------------------------------------------------- 660 661 IF( ln_abl ) THEN !== ABL formulation ==! multiplication by rho_air and turbulent fluxes computation done in ablstp 662 !! FL do we need this multiplication by tmask ... ??? 663 DO jj = 1, jpj 664 DO ji = 1, jpi 665 zztmp = zU_zu(ji,jj) !* tmask(ji,jj,1) 666 wndm(ji,jj) = zztmp ! Store zU_zu in wndm to compute ustar2 in ablmod 667 pcd_du(ji,jj) = zztmp * zcd_oce(ji,jj) 668 psen(ji,jj) = zztmp * zch_oce(ji,jj) 669 pevp(ji,jj) = zztmp * zce_oce(ji,jj) 670 END DO 671 END DO 672 ELSE !== BLK formulation ==! turbulent fluxes computation 673 CALL BULK_FORMULA( rn_zu, zst(:,:), pssq(:,:), t_zu(:,:), q_zu(:,:), & 674 & zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:), & 675 & wndm(:,:), zU_zu(:,:), pslp(:,:), & 676 & taum(:,:), psen(:,:), zqla(:,:), & 677 & pEvap=pevp(:,:), prhoa=rhoa(:,:) ) 678 679 zqla(:,:) = zqla(:,:) * tmask(:,:,1) 680 psen(:,:) = psen(:,:) * tmask(:,:,1) 681 taum(:,:) = taum(:,:) * tmask(:,:,1) 682 pevp(:,:) = pevp(:,:) * tmask(:,:,1) 683 684 ! Tau i and j component on T-grid points, using array "zcd_oce" as a temporary array... 685 zcd_oce = 0._wp 686 WHERE ( wndm > 0._wp ) zcd_oce = taum / wndm 687 zwnd_i = zcd_oce * zwnd_i 688 zwnd_j = zcd_oce * zwnd_j 689 690 CALL iom_put( "taum_oce", taum ) ! output wind stress module 691 692 ! ... utau, vtau at U- and V_points, resp. 693 ! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines 694 ! Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves 695 DO jj = 1, jpjm1 696 DO ji = 1, fs_jpim1 697 utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj ) ) & 698 & * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1)) 699 vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji ,jj+1) ) & 700 & * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1)) 701 END DO 702 END DO 703 CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. ) 704 705 IF(ln_ctl) THEN 706 CALL prt_ctl( tab2d_1=wndm , clinfo1=' blk_oce_1: wndm : ') 707 CALL prt_ctl( tab2d_1=utau , clinfo1=' blk_oce_1: utau : ', mask1=umask, & 708 & tab2d_2=vtau , clinfo2=' vtau : ', mask2=vmask ) 709 ENDIF 710 ! 711 ENDIF 712 ! 713 IF(ln_ctl) THEN 714 CALL prt_ctl( tab2d_1=pevp , clinfo1=' blk_oce_1: pevp : ' ) 715 CALL prt_ctl( tab2d_1=psen , clinfo1=' blk_oce_1: psen : ' ) 716 CALL prt_ctl( tab2d_1=pssq , clinfo1=' blk_oce_1: pssq : ' ) 717 ENDIF 718 ! 719 END SUBROUTINE blk_oce_1 720 721 722 SUBROUTINE blk_oce_2( ptair, pqsr, pqlw, pprec, & ! <<= in 723 & psnow, pst , psen, pevp ) ! <<= in 724 !!--------------------------------------------------------------------- 725 !! *** ROUTINE blk_oce_2 *** 726 !! 727 !! ** Purpose : finalize the momentum, heat and freshwater fluxes computation 728 !! at the ocean surface at each time step knowing Cd, Ch, Ce and 729 !! atmospheric variables (from ABL or external data) 354 730 !! 355 731 !! ** Outputs : - utau : i-component of the stress at U-point (N/m2) … … 360 736 !! - qns : Non Solar heat flux over the ocean (W/m2) 361 737 !! - emp : evaporation minus precipitation (kg/m2/s) 362 !! 363 !! ** Nota : sf has to be a dummy argument for AGRIF on NEC 364 !!--------------------------------------------------------------------- 365 INTEGER , INTENT(in ) :: kt ! time step index 366 TYPE(fld), INTENT(inout), DIMENSION(:) :: sf ! input data 367 REAL(wp) , INTENT(in) , DIMENSION(:,:) :: pst ! surface temperature [Celcius] 368 REAL(wp) , INTENT(in) , DIMENSION(:,:) :: pu ! surface current at U-point (i-component) [m/s] 369 REAL(wp) , INTENT(in) , DIMENSION(:,:) :: pv ! surface current at V-point (j-component) [m/s] 738 !!--------------------------------------------------------------------- 739 REAL(wp), INTENT(in), DIMENSION(:,:) :: ptair 740 REAL(wp), INTENT(in), DIMENSION(:,:) :: pqsr 741 REAL(wp), INTENT(in), DIMENSION(:,:) :: pqlw 742 REAL(wp), INTENT(in), DIMENSION(:,:) :: pprec 743 REAL(wp), INTENT(in), DIMENSION(:,:) :: psnow 744 REAL(wp), INTENT(in), DIMENSION(:,:) :: pst ! surface temperature [Celcius] 745 REAL(wp), INTENT(in), DIMENSION(:,:) :: psen 746 REAL(wp), INTENT(in), DIMENSION(:,:) :: pevp 370 747 ! 371 748 INTEGER :: ji, jj ! dummy loop indices 372 REAL(wp) :: zztmp ! local variable 373 REAL(wp), DIMENSION(jpi,jpj) :: zwnd_i, zwnd_j ! wind speed components at T-point 374 REAL(wp), DIMENSION(jpi,jpj) :: zsq ! specific humidity at pst 375 REAL(wp), DIMENSION(jpi,jpj) :: zqlw, zqsb ! long wave and sensible heat fluxes 376 REAL(wp), DIMENSION(jpi,jpj) :: zqla, zevap ! latent heat fluxes and evaporation 749 REAL(wp) :: zztmp,zz1,zz2,zz3 ! local variable 750 REAL(wp), DIMENSION(jpi,jpj) :: zqlw ! long wave and sensible heat fluxes 751 REAL(wp), DIMENSION(jpi,jpj) :: zqla ! latent heat fluxes and evaporation 377 752 REAL(wp), DIMENSION(jpi,jpj) :: zst ! surface temperature in Kelvin 378 REAL(wp), DIMENSION(jpi,jpj) :: zU_zu ! bulk wind speed at height zu [m/s]379 REAL(wp), DIMENSION(jpi,jpj) :: ztpot ! potential temperature of air at z=rn_zqt [K]380 REAL(wp), DIMENSION(jpi,jpj) :: zrhoa ! density of air [kg/m^3]381 753 !!--------------------------------------------------------------------- 382 754 ! … … 384 756 zst(:,:) = pst(:,:) + rt0 ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K) 385 757 758 386 759 ! ----------------------------------------------------------------------------- ! 387 ! 0 Wind components and module at T-point relative to the moving ocean!760 ! III Net longwave radiative FLUX ! 388 761 ! ----------------------------------------------------------------------------- ! 389 762 390 ! ... components ( U10m - U_oce ) at T-point (unmasked) 391 !!gm move zwnd_i (_j) set to zero inside the key_cyclone ??? 392 zwnd_i(:,:) = 0._wp 393 zwnd_j(:,:) = 0._wp 394 #if defined key_cyclone 395 CALL wnd_cyc( kt, zwnd_i, zwnd_j ) ! add analytical tropical cyclone (Vincent et al. JGR 2012) 396 DO jj = 2, jpjm1 397 DO ji = fs_2, fs_jpim1 ! vect. opt. 398 sf(jp_wndi)%fnow(ji,jj,1) = sf(jp_wndi)%fnow(ji,jj,1) + zwnd_i(ji,jj) 399 sf(jp_wndj)%fnow(ji,jj,1) = sf(jp_wndj)%fnow(ji,jj,1) + zwnd_j(ji,jj) 400 END DO 401 END DO 402 #endif 403 DO jj = 2, jpjm1 404 DO ji = fs_2, fs_jpim1 ! vect. opt. 405 zwnd_i(ji,jj) = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj ) + pu(ji,jj) ) ) 406 zwnd_j(ji,jj) = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji ,jj-1) + pv(ji,jj) ) ) 407 END DO 408 END DO 409 CALL lbc_lnk_multi( 'sbcblk', zwnd_i, 'T', -1., zwnd_j, 'T', -1. ) 410 ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked) 411 wndm(:,:) = SQRT( zwnd_i(:,:) * zwnd_i(:,:) & 412 & + zwnd_j(:,:) * zwnd_j(:,:) ) * tmask(:,:,1) 413 414 ! ----------------------------------------------------------------------------- ! 415 ! I Radiative FLUXES ! 416 ! ----------------------------------------------------------------------------- ! 417 418 ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle ! Short Wave 419 zztmp = 1. - albo 420 IF( ln_dm2dc ) THEN ; qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1) 421 ELSE ; qsr(:,:) = zztmp * sf(jp_qsr)%fnow(:,:,1) * tmask(:,:,1) 422 ENDIF 423 424 zqlw(:,:) = ( sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:) ) * tmask(:,:,1) ! Long Wave 425 426 ! ----------------------------------------------------------------------------- ! 427 ! II Turbulent FLUXES ! 428 ! ----------------------------------------------------------------------------- ! 429 430 ! ... specific humidity at SST and IST tmask( 431 zsq(:,:) = 0.98 * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) ) 432 !! 433 !! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate 434 !! (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2 435 !! (since reanalysis products provide T at z, not theta !) 436 ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1) ) * rn_zqt 437 438 SELECT CASE( nblk ) !== transfer coefficients ==! Cd, Ch, Ce at T-point 439 ! 440 CASE( np_NCAR ) ; CALL turb_ncar ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! NCAR-COREv2 441 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 442 CASE( np_COARE_3p0 ) ; CALL turb_coare ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! COARE v3.0 443 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 444 CASE( np_COARE_3p5 ) ; CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! COARE v3.5 445 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 446 CASE( np_ECMWF ) ; CALL turb_ecmwf ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm, & ! ECMWF 447 & Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 448 CASE DEFAULT 449 CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 450 END SELECT 451 452 ! ! Compute true air density : 453 IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN ! At zu: (probably useless to remove zrho*grav*rn_zu from SLP...) 454 zrhoa(:,:) = rho_air( t_zu(:,:) , q_zu(:,:) , sf(jp_slp)%fnow(:,:,1) ) 455 ELSE ! At zt: 456 zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 457 END IF 458 459 !! CALL iom_put( "Cd_oce", Cd_atm) ! output value of pure ocean-atm. transfer coef. 460 !! CALL iom_put( "Ch_oce", Ch_atm) ! output value of pure ocean-atm. transfer coef. 461 462 DO jj = 1, jpj ! tau module, i and j component 463 DO ji = 1, jpi 464 zztmp = zrhoa(ji,jj) * zU_zu(ji,jj) * Cd_atm(ji,jj) ! using bulk wind speed 465 taum (ji,jj) = zztmp * wndm (ji,jj) 466 zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 467 zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 468 END DO 469 END DO 470 471 ! ! add the HF tau contribution to the wind stress module 472 IF( lhftau ) taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1) 473 474 CALL iom_put( "taum_oce", taum ) ! output wind stress module 475 476 ! ... utau, vtau at U- and V_points, resp. 477 ! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines 478 ! Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves 479 DO jj = 1, jpjm1 480 DO ji = 1, fs_jpim1 481 utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj ) ) & 482 & * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1)) 483 vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji ,jj+1) ) & 484 & * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1)) 485 END DO 486 END DO 487 CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. ) 763 !! LB: now moved after Turbulent fluxes because must use the skin temperature rather that the SST 764 !! (zst is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.) 765 zqlw(:,:) = emiss_w * ( pqlw(:,:) - stefan*zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:) ) * tmask(:,:,1) ! Net radiative longwave flux 488 766 489 767 ! Turbulent fluxes over ocean 490 768 ! ----------------------------- 491 769 492 ! zqla used as temporary array, for rho*U (common term of bulk formulae): 493 zqla(:,:) = zrhoa(:,:) * zU_zu(:,:) * tmask(:,:,1) 494 495 IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN 496 !! q_air and t_air are given at 10m (wind reference height) 497 zevap(:,:) = rn_efac*MAX( 0._wp, zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed 498 zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:) ) ! Sensible Heat, using bulk wind speed 499 ELSE 500 !! q_air and t_air are not given at 10m (wind reference height) 501 ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!! 502 zevap(:,:) = rn_efac*MAX( 0._wp, zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - q_zu(:,:) ) ) ! Evaporation, using bulk wind speed 503 zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) ) ! Sensible Heat, using bulk wind speed 504 ENDIF 505 506 zqla(:,:) = L_vap(zst(:,:)) * zevap(:,:) ! Latent Heat flux 507 770 ! use scalar version of L_vap() for AGRIF compatibility 771 DO jj = 1, jpj 772 DO ji = 1, jpi 773 zqla(ji,jj) = L_vap( zst(ji,jj) ) * pevp(ji,jj) * -1._wp ! Latent Heat flux !!GS: possibility to add a global qla to avoid recomputation after abl update 774 ENDDO 775 ENDDO 508 776 509 777 IF(ln_ctl) THEN 510 CALL prt_ctl( tab2d_1=zqla , clinfo1=' blk_oce: zqla : ', tab2d_2=Ce_atm , clinfo2=' Ce_oce : ' ) 511 CALL prt_ctl( tab2d_1=zqsb , clinfo1=' blk_oce: zqsb : ', tab2d_2=Ch_atm , clinfo2=' Ch_oce : ' ) 512 CALL prt_ctl( tab2d_1=zqlw , clinfo1=' blk_oce: zqlw : ', tab2d_2=qsr, clinfo2=' qsr : ' ) 513 CALL prt_ctl( tab2d_1=zsq , clinfo1=' blk_oce: zsq : ', tab2d_2=zst, clinfo2=' zst : ' ) 514 CALL prt_ctl( tab2d_1=utau , clinfo1=' blk_oce: utau : ', mask1=umask, & 515 & tab2d_2=vtau , clinfo2= ' vtau : ', mask2=vmask ) 516 CALL prt_ctl( tab2d_1=wndm , clinfo1=' blk_oce: wndm : ') 517 CALL prt_ctl( tab2d_1=zst , clinfo1=' blk_oce: zst : ') 778 CALL prt_ctl( tab2d_1=zqla , clinfo1=' blk_oce_2: zqla : ' ) 779 CALL prt_ctl( tab2d_1=zqlw , clinfo1=' blk_oce_2: zqlw : ', tab2d_2=qsr, clinfo2=' qsr : ' ) 780 518 781 ENDIF 519 782 520 783 ! ----------------------------------------------------------------------------- ! 521 ! I IITotal FLUXES !784 ! IV Total FLUXES ! 522 785 ! ----------------------------------------------------------------------------- ! 523 786 ! 524 emp (:,:) = ( zevap(:,:)& ! mass flux (evap. - precip.)525 & - sf(jp_prec)%fnow(:,:,1) * rn_pfac ) * tmask(:,:,1)526 ! 527 qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)& ! Downward Non Solar528 & - sf(jp_snow)%fnow(:,:,1) * rn_pfac * rLfus & ! remove latent melting heat for solid precip529 & - zevap(:,:) * pst(:,:) * rcp & ! remove evap heat content at SST530 & + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) ) * rn_pfac& ! add liquid precip heat content at Tair531 & * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp &532 & + sf(jp_snow)%fnow(:,:,1) * rn_pfac & ! add solid precip heat content at min(Tair,Tsnow)533 & * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi787 emp (:,:) = ( pevp(:,:) & ! mass flux (evap. - precip.) 788 & - pprec(:,:) * rn_pfac ) * tmask(:,:,1) 789 ! 790 qns(:,:) = zqlw(:,:) + psen(:,:) + zqla(:,:) & ! Downward Non Solar 791 & - psnow(:,:) * rn_pfac * rLfus & ! remove latent melting heat for solid precip 792 & - pevp(:,:) * pst(:,:) * rcp & ! remove evap heat content at SST !LB??? pst is Celsius !? 793 & + ( pprec(:,:) - psnow(:,:) ) * rn_pfac & ! add liquid precip heat content at Tair 794 & * ( ptair(:,:) - rt0 ) * rcp & 795 & + psnow(:,:) * rn_pfac & ! add solid precip heat content at min(Tair,Tsnow) 796 & * ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi 534 797 qns(:,:) = qns(:,:) * tmask(:,:,1) 535 798 ! 536 799 #if defined key_si3 537 qns_oce(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)! non solar without emp (only needed by SI3)800 qns_oce(:,:) = zqlw(:,:) + psen(:,:) + zqla(:,:) ! non solar without emp (only needed by SI3) 538 801 qsr_oce(:,:) = qsr(:,:) 539 802 #endif 540 803 ! 804 CALL iom_put( "rho_air" , rhoa*tmask(:,:,1) ) ! output air density [kg/m^3] 805 CALL iom_put( "evap_oce" , pevp ) ! evaporation 806 CALL iom_put( "qlw_oce" , zqlw ) ! output downward longwave heat over the ocean 807 CALL iom_put( "qsb_oce" , psen ) ! output downward sensible heat over the ocean 808 CALL iom_put( "qla_oce" , zqla ) ! output downward latent heat over the ocean 809 tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s] 810 sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s] 811 CALL iom_put( 'snowpre', sprecip ) ! Snow 812 CALL iom_put( 'precip' , tprecip ) ! Total precipitation 813 ! 541 814 IF ( nn_ice == 0 ) THEN 542 CALL iom_put( "qlw_oce" , zqlw ) ! output downward longwave heat over the ocean 543 CALL iom_put( "qsb_oce" , - zqsb ) ! output downward sensible heat over the ocean 544 CALL iom_put( "qla_oce" , - zqla ) ! output downward latent heat over the ocean 545 CALL iom_put( "qemp_oce", qns-zqlw+zqsb+zqla ) ! output downward heat content of E-P over the ocean 546 CALL iom_put( "qns_oce" , qns ) ! output downward non solar heat over the ocean 547 CALL iom_put( "qsr_oce" , qsr ) ! output downward solar heat over the ocean 548 CALL iom_put( "qt_oce" , qns+qsr ) ! output total downward heat over the ocean 549 tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s] 550 sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s] 551 CALL iom_put( 'snowpre', sprecip ) ! Snow 552 CALL iom_put( 'precip' , tprecip ) ! Total precipitation 815 CALL iom_put( "qemp_oce" , qns-zqlw-psen-zqla ) ! output downward heat content of E-P over the ocean 816 CALL iom_put( "qns_oce" , qns ) ! output downward non solar heat over the ocean 817 CALL iom_put( "qsr_oce" , qsr ) ! output downward solar heat over the ocean 818 CALL iom_put( "qt_oce" , qns+qsr ) ! output total downward heat over the ocean 819 ENDIF 820 ! 821 IF( ln_skin_cs .OR. ln_skin_wl ) THEN 822 CALL iom_put( "t_skin" , (zst - rt0) * tmask(:,:,1) ) ! T_skin in Celsius 823 CALL iom_put( "dt_skin" , (zst - pst - rt0) * tmask(:,:,1) ) ! T_skin - SST temperature difference... 553 824 ENDIF 554 825 ! 555 826 IF(ln_ctl) THEN 556 CALL prt_ctl(tab2d_1=zqsb , clinfo1=' blk_oce: zqsb : ', tab2d_2=zqlw , clinfo2=' zqlw : ') 557 CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce: zqla : ', tab2d_2=qsr , clinfo2=' qsr : ') 558 CALL prt_ctl(tab2d_1=pst , clinfo1=' blk_oce: pst : ', tab2d_2=emp , clinfo2=' emp : ') 559 CALL prt_ctl(tab2d_1=utau , clinfo1=' blk_oce: utau : ', mask1=umask, & 560 & tab2d_2=vtau , clinfo2= ' vtau : ' , mask2=vmask ) 561 ENDIF 562 ! 563 END SUBROUTINE blk_oce 564 565 566 567 FUNCTION rho_air( ptak, pqa, pslp ) 568 !!------------------------------------------------------------------------------- 569 !! *** FUNCTION rho_air *** 570 !! 571 !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere 572 !! 573 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 574 !!------------------------------------------------------------------------------- 575 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K] 576 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg] 577 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp ! pressure in [Pa] 578 REAL(wp), DIMENSION(jpi,jpj) :: rho_air ! density of moist air [kg/m^3] 579 !!------------------------------------------------------------------------------- 580 ! 581 rho_air = pslp / ( R_dry*ptak * ( 1._wp + rctv0*pqa ) ) 582 ! 583 END FUNCTION rho_air 584 585 586 FUNCTION cp_air( pqa ) 587 !!------------------------------------------------------------------------------- 588 !! *** FUNCTION cp_air *** 589 !! 590 !! ** Purpose : Compute specific heat (Cp) of moist air 591 !! 592 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 593 !!------------------------------------------------------------------------------- 594 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg] 595 REAL(wp), DIMENSION(jpi,jpj) :: cp_air ! specific heat of moist air [J/K/kg] 596 !!------------------------------------------------------------------------------- 597 ! 598 Cp_air = Cp_dry + Cp_vap * pqa 599 ! 600 END FUNCTION cp_air 601 602 603 FUNCTION q_sat( ptak, pslp ) 604 !!---------------------------------------------------------------------------------- 605 !! *** FUNCTION q_sat *** 606 !! 607 !! ** Purpose : Specific humidity at saturation in [kg/kg] 608 !! Based on accurate estimate of "e_sat" 609 !! aka saturation water vapor (Goff, 1957) 610 !! 611 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 612 !!---------------------------------------------------------------------------------- 613 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K] 614 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp ! sea level atmospheric pressure [Pa] 615 REAL(wp), DIMENSION(jpi,jpj) :: q_sat ! Specific humidity at saturation [kg/kg] 616 ! 617 INTEGER :: ji, jj ! dummy loop indices 618 REAL(wp) :: ze_sat, ztmp ! local scalar 619 !!---------------------------------------------------------------------------------- 620 ! 621 DO jj = 1, jpj 622 DO ji = 1, jpi 623 ! 624 ztmp = rt0 / ptak(ji,jj) 625 ! 626 ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957) 627 ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0) & 628 & + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) ) & 629 & + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614 ) 630 ! 631 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] 632 ! 633 END DO 634 END DO 635 ! 636 END FUNCTION q_sat 637 638 639 FUNCTION gamma_moist( ptak, pqa ) 640 !!---------------------------------------------------------------------------------- 641 !! *** FUNCTION gamma_moist *** 642 !! 643 !! ** Purpose : Compute the moist adiabatic lapse-rate. 644 !! => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate 645 !! => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html 646 !! 647 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 648 !!---------------------------------------------------------------------------------- 649 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K] 650 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! specific humidity [kg/kg] 651 REAL(wp), DIMENSION(jpi,jpj) :: gamma_moist ! moist adiabatic lapse-rate 652 ! 653 INTEGER :: ji, jj ! dummy loop indices 654 REAL(wp) :: zrv, ziRT ! local scalar 655 !!---------------------------------------------------------------------------------- 656 ! 657 DO jj = 1, jpj 658 DO ji = 1, jpi 659 zrv = pqa(ji,jj) / (1. - pqa(ji,jj)) 660 ziRT = 1. / (R_dry*ptak(ji,jj)) ! 1/RT 661 gamma_moist(ji,jj) = grav * ( 1. + rLevap*zrv*ziRT ) / ( Cp_dry + rLevap*rLevap*zrv*reps0*ziRT/ptak(ji,jj) ) 662 END DO 663 END DO 664 ! 665 END FUNCTION gamma_moist 666 667 668 FUNCTION L_vap( psst ) 669 !!--------------------------------------------------------------------------------- 670 !! *** FUNCTION L_vap *** 671 !! 672 !! ** Purpose : Compute the latent heat of vaporization of water from temperature 673 !! 674 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 675 !!---------------------------------------------------------------------------------- 676 REAL(wp), DIMENSION(jpi,jpj) :: L_vap ! latent heat of vaporization [J/kg] 677 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst ! water temperature [K] 678 !!---------------------------------------------------------------------------------- 679 ! 680 L_vap = ( 2.501 - 0.00237 * ( psst(:,:) - rt0) ) * 1.e6 681 ! 682 END FUNCTION L_vap 827 CALL prt_ctl(tab2d_1=zqlw , clinfo1=' blk_oce_2: zqlw : ') 828 CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce_2: zqla : ', tab2d_2=qsr , clinfo2=' qsr : ') 829 CALL prt_ctl(tab2d_1=emp , clinfo1=' blk_oce_2: emp : ') 830 ENDIF 831 ! 832 END SUBROUTINE blk_oce_2 833 683 834 684 835 #if defined key_si3 … … 686 837 !! 'key_si3' SI3 sea-ice model 687 838 !!---------------------------------------------------------------------- 688 !! blk_ice_ tau: provide the air-ice stress689 !! blk_ice_ flx: provide the heat and mass fluxes at air-ice interface839 !! blk_ice_1 : provide the air-ice stress 840 !! blk_ice_2 : provide the heat and mass fluxes at air-ice interface 690 841 !! blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux) 691 842 !! Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag 692 !! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 843 !! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 693 844 !!---------------------------------------------------------------------- 694 845 695 SUBROUTINE blk_ice_tau 696 !!--------------------------------------------------------------------- 697 !! *** ROUTINE blk_ice_tau *** 846 SUBROUTINE blk_ice_1( pwndi, pwndj, ptair, phumi, pslp , puice, pvice, ptsui, & ! inputs 847 & putaui, pvtaui, pseni, pevpi, pssqi, pcd_dui ) ! optional outputs 848 !!--------------------------------------------------------------------- 849 !! *** ROUTINE blk_ice_1 *** 698 850 !! 699 851 !! ** Purpose : provide the surface boundary condition over sea-ice … … 703 855 !! NB: ice drag coefficient is assumed to be a constant 704 856 !!--------------------------------------------------------------------- 857 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pslp ! sea-level pressure [Pa] 858 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndi ! atmospheric wind at T-point [m/s] 859 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pwndj ! atmospheric wind at T-point [m/s] 860 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: ptair ! atmospheric wind at T-point [m/s] 861 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: phumi ! atmospheric wind at T-point [m/s] 862 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: puice ! sea-ice velocity on I or C grid [m/s] 863 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: pvice ! " 864 REAL(wp) , INTENT(in ), DIMENSION(:,: ) :: ptsui ! sea-ice surface temperature [K] 865 REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: putaui ! if ln_blk 866 REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pvtaui ! if ln_blk 867 REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pseni ! if ln_abl 868 REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pevpi ! if ln_abl 869 REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pssqi ! if ln_abl 870 REAL(wp) , INTENT( out), DIMENSION(:,: ), OPTIONAL :: pcd_dui ! if ln_abl 871 ! 705 872 INTEGER :: ji, jj ! dummy loop indices 706 REAL(wp) :: zwndi_f , zwndj_f, zwnorm_f ! relative wind module and components at F-point707 873 REAL(wp) :: zwndi_t , zwndj_t ! relative wind components at T-point 708 REAL(wp), DIMENSION(jpi,jpj) :: zrhoa ! transfer coefficient for momentum (tau) 709 !!--------------------------------------------------------------------- 710 ! 711 ! set transfer coefficients to default sea-ice values 712 Cd_atm(:,:) = Cd_ice 713 Ch_atm(:,:) = Cd_ice 714 Ce_atm(:,:) = Cd_ice 715 716 wndm_ice(:,:) = 0._wp !!gm brutal.... 874 REAL(wp) :: zootm_su ! sea-ice surface mean temperature 875 REAL(wp) :: zztmp1, zztmp2 ! temporary arrays 876 REAL(wp), DIMENSION(jpi,jpj) :: zcd_dui ! transfer coefficient for momentum (tau) 877 !!--------------------------------------------------------------------- 878 ! 717 879 718 880 ! ------------------------------------------------------------ ! … … 722 884 DO jj = 2, jpjm1 723 885 DO ji = fs_2, fs_jpim1 ! vect. opt. 724 zwndi_t = ( sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) )725 zwndj_t = ( sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji ,jj-1) + v_ice(ji,jj) ) )886 zwndi_t = ( pwndi(ji,jj) - rn_vfac * 0.5_wp * ( puice(ji-1,jj ) + puice(ji,jj) ) ) 887 zwndj_t = ( pwndj(ji,jj) - rn_vfac * 0.5_wp * ( pvice(ji ,jj-1) + pvice(ji,jj) ) ) 726 888 wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 727 889 END DO … … 731 893 ! Make ice-atm. drag dependent on ice concentration 732 894 IF ( ln_Cd_L12 ) THEN ! calculate new drag from Lupkes(2012) equations 733 CALL Cdn10_Lupkes2012( Cd_atm ) 734 Ch_atm(:,:) = Cd_atm(:,:) ! momentum and heat transfer coef. are considered identical 895 CALL Cdn10_Lupkes2012( Cd_ice ) 896 Ch_ice(:,:) = Cd_ice(:,:) ! momentum and heat transfer coef. are considered identical 897 Ce_ice(:,:) = Cd_ice(:,:) 735 898 ELSEIF( ln_Cd_L15 ) THEN ! calculate new drag from Lupkes(2015) equations 736 CALL Cdn10_Lupkes2015( Cd_atm, Ch_atm ) 737 ENDIF 738 739 !! CALL iom_put( "Cd_ice", Cd_atm) ! output value of pure ice-atm. transfer coef. 740 !! CALL iom_put( "Ch_ice", Ch_atm) ! output value of pure ice-atm. transfer coef. 899 CALL Cdn10_Lupkes2015( ptsui, pslp, Cd_ice, Ch_ice ) 900 Ce_ice(:,:) = Ch_ice(:,:) ! sensible and latent heat transfer coef. are considered identical 901 ENDIF 902 903 !! IF ( iom_use("Cd_ice") ) CALL iom_put("Cd_ice", Cd_ice) ! output value of pure ice-atm. transfer coef. 904 !! IF ( iom_use("Ch_ice") ) CALL iom_put("Ch_ice", Ch_ice) ! output value of pure ice-atm. transfer coef. 741 905 742 906 ! local scalars ( place there for vector optimisation purposes) 743 ! Computing density of air! Way denser that 1.2 over sea-ice !!! 744 zrhoa (:,:) = rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) 745 746 !!gm brutal.... 747 utau_ice (:,:) = 0._wp 748 vtau_ice (:,:) = 0._wp 749 !!gm end 750 751 ! ------------------------------------------------------------ ! 752 ! Wind stress relative to the moving ice ( U10m - U_ice ) ! 753 ! ------------------------------------------------------------ ! 754 ! C-grid ice dynamics : U & V-points (same as ocean) 755 DO jj = 2, jpjm1 756 DO ji = fs_2, fs_jpim1 ! vect. opt. 757 utau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji+1,jj ) + wndm_ice(ji,jj) ) & 758 & * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) ) 759 vtau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji,jj+1 ) + wndm_ice(ji,jj) ) & 760 & * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) ) 907 !IF (ln_abl) rhoa (:,:) = rho_air( ptair(:,:), phumi(:,:), pslp(:,:) ) !!GS: rhoa must be (re)computed here with ABL to avoid division by zero after (TBI) 908 zcd_dui(:,:) = wndm_ice(:,:) * Cd_ice(:,:) 909 910 IF( ln_blk ) THEN 911 ! ------------------------------------------------------------ ! 912 ! Wind stress relative to the moving ice ( U10m - U_ice ) ! 913 ! ------------------------------------------------------------ ! 914 ! C-grid ice dynamics : U & V-points (same as ocean) 915 DO jj = 2, jpjm1 916 DO ji = fs_2, fs_jpim1 ! vect. opt. 917 putaui(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * zcd_dui(ji+1,jj) & 918 & + rhoa(ji ,jj) * zcd_dui(ji ,jj) ) & 919 & * ( 0.5_wp * ( pwndi(ji+1,jj) + pwndi(ji,jj) ) - rn_vfac * puice(ji,jj) ) 920 pvtaui(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * zcd_dui(ji,jj+1) & 921 & + rhoa(ji,jj ) * zcd_dui(ji,jj ) ) & 922 & * ( 0.5_wp * ( pwndj(ji,jj+1) + pwndj(ji,jj) ) - rn_vfac * pvice(ji,jj) ) 923 END DO 761 924 END DO 762 END DO 763 CALL lbc_lnk_multi( 'sbcblk', utau_ice, 'U', -1., vtau_ice, 'V', -1. ) 764 ! 765 ! 766 IF(ln_ctl) THEN 767 CALL prt_ctl(tab2d_1=utau_ice , clinfo1=' blk_ice: utau_ice : ', tab2d_2=vtau_ice , clinfo2=' vtau_ice : ') 768 CALL prt_ctl(tab2d_1=wndm_ice , clinfo1=' blk_ice: wndm_ice : ') 769 ENDIF 770 ! 771 END SUBROUTINE blk_ice_tau 772 773 774 SUBROUTINE blk_ice_flx( ptsu, phs, phi, palb ) 775 !!--------------------------------------------------------------------- 776 !! *** ROUTINE blk_ice_flx *** 925 CALL lbc_lnk_multi( 'sbcblk', putaui, 'U', -1., pvtaui, 'V', -1. ) 926 ! 927 IF(ln_ctl) CALL prt_ctl( tab2d_1=putaui , clinfo1=' blk_ice: putaui : ' & 928 & , tab2d_2=pvtaui , clinfo2=' pvtaui : ' ) 929 ELSE 930 zztmp1 = 11637800.0_wp 931 zztmp2 = -5897.8_wp 932 DO jj = 1, jpj 933 DO ji = 1, jpi 934 pcd_dui(ji,jj) = zcd_dui (ji,jj) 935 pseni (ji,jj) = wndm_ice(ji,jj) * Ch_ice(ji,jj) 936 pevpi (ji,jj) = wndm_ice(ji,jj) * Ce_ice(ji,jj) 937 zootm_su = zztmp2 / ptsui(ji,jj) ! ptsui is in K (it can't be zero ??) 938 pssqi (ji,jj) = zztmp1 * EXP( zootm_su ) / rhoa(ji,jj) 939 END DO 940 END DO 941 ENDIF 942 ! 943 IF(ln_ctl) CALL prt_ctl(tab2d_1=wndm_ice , clinfo1=' blk_ice: wndm_ice : ') 944 ! 945 END SUBROUTINE blk_ice_1 946 947 948 SUBROUTINE blk_ice_2( ptsu, phs, phi, palb, ptair, phumi, pslp, pqlw, pprec, psnow ) 949 !!--------------------------------------------------------------------- 950 !! *** ROUTINE blk_ice_2 *** 777 951 !! 778 952 !! ** Purpose : provide the heat and mass fluxes at air-ice interface … … 784 958 !! caution : the net upward water flux has with mm/day unit 785 959 !!--------------------------------------------------------------------- 786 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: ptsu ! sea ice surface temperature 960 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: ptsu ! sea ice surface temperature [K] 787 961 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phs ! snow thickness 788 962 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: phi ! ice thickness 789 963 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: palb ! ice albedo (all skies) 964 REAL(wp), DIMENSION(:,: ), INTENT(in) :: ptair 965 REAL(wp), DIMENSION(:,: ), INTENT(in) :: phumi 966 REAL(wp), DIMENSION(:,: ), INTENT(in) :: pslp 967 REAL(wp), DIMENSION(:,: ), INTENT(in) :: pqlw 968 REAL(wp), DIMENSION(:,: ), INTENT(in) :: pprec 969 REAL(wp), DIMENSION(:,: ), INTENT(in) :: psnow 790 970 !! 791 971 INTEGER :: ji, jj, jl ! dummy loop indices 792 972 REAL(wp) :: zst3 ! local variable 793 973 REAL(wp) :: zcoef_dqlw, zcoef_dqla ! - - 794 REAL(wp) :: zztmp, z 1_rLsub! - -974 REAL(wp) :: zztmp, zztmp2, z1_rLsub ! - - 795 975 REAL(wp) :: zfr1, zfr2 ! local variables 796 976 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_st ! inverse of surface temperature … … 800 980 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqsb ! sensible heat sensitivity over ice 801 981 REAL(wp), DIMENSION(jpi,jpj) :: zevap, zsnw ! evaporation and snw distribution after wind blowing (SI3) 802 REAL(wp), DIMENSION(jpi,jpj) :: z rhoa982 REAL(wp), DIMENSION(jpi,jpj) :: zqair ! specific humidity of air at z=rn_zqt [kg/kg] !LB 803 983 REAL(wp), DIMENSION(jpi,jpj) :: ztmp, ztmp2 804 984 !!--------------------------------------------------------------------- 805 985 ! 806 zcoef_dqlw = 4.0 * 0.95 * Stef ! local scalars 807 zcoef_dqla = -Ls * 11637800. * (-5897.8) 808 ! 809 zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 986 zcoef_dqlw = 4._wp * 0.95_wp * stefan ! local scalars 987 zcoef_dqla = -rLsub * 11637800._wp * (-5897.8_wp) !LB: BAD! 988 ! 989 SELECT CASE( nhumi ) 990 CASE( np_humi_sph ) 991 zqair(:,:) = phumi(:,:) ! what we read in file is already a spec. humidity! 992 CASE( np_humi_dpt ) 993 zqair(:,:) = q_sat( phumi(:,:), pslp ) 994 CASE( np_humi_rlh ) 995 zqair(:,:) = q_air_rh( 0.01_wp*phumi(:,:), ptair(:,:), pslp(:,:) ) !LB: 0.01 => RH is % percent in file 996 END SELECT 810 997 ! 811 998 zztmp = 1. / ( 1. - albo ) 812 WHERE( ptsu(:,:,:) /= 0._wp ) ; z1_st(:,:,:) = 1._wp / ptsu(:,:,:) 813 ELSEWHERE ; z1_st(:,:,:) = 0._wp 999 WHERE( ptsu(:,:,:) /= 0._wp ) 1000 z1_st(:,:,:) = 1._wp / ptsu(:,:,:) 1001 ELSEWHERE 1002 z1_st(:,:,:) = 0._wp 814 1003 END WHERE 815 1004 ! ! ========================== ! … … 825 1014 qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj) 826 1015 ! Long Wave (lw) 827 z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef* ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)1016 z_qlw(ji,jj,jl) = 0.95 * ( pqlw(ji,jj) - stefan * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 828 1017 ! lw sensitivity 829 1018 z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3 … … 833 1022 ! ----------------------------! 834 1023 835 ! ... turbulent heat fluxes with Ch_ atm recalculated in blk_ice_tau1024 ! ... turbulent heat fluxes with Ch_ice recalculated in blk_ice_1 836 1025 ! Sensible Heat 837 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))1026 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)) 838 1027 ! Latent Heat 839 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls * Ch_atm(ji,jj) * wndm_ice(ji,jj) * & 840 & ( 11637800. * EXP( -5897.8 * z1_st(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) ) 1028 zztmp2 = EXP( -5897.8 * z1_st(ji,jj,jl) ) 1029 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa(ji,jj) * rLsub * Ce_ice(ji,jj) * wndm_ice(ji,jj) * & 1030 & ( 11637800. * zztmp2 / rhoa(ji,jj) - zqair(ji,jj) ) ) 841 1031 ! Latent heat sensitivity for ice (Dqla/Dt) 842 1032 IF( qla_ice(ji,jj,jl) > 0._wp ) THEN 843 dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * C h_atm(ji,jj) * wndm_ice(ji,jj) * &844 & z1_st(ji,jj,jl) *z1_st(ji,jj,jl) * EXP(-5897.8 * z1_st(ji,jj,jl))1033 dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ce_ice(ji,jj) * wndm_ice(ji,jj) * & 1034 & z1_st(ji,jj,jl) * z1_st(ji,jj,jl) * zztmp2 845 1035 ELSE 846 1036 dqla_ice(ji,jj,jl) = 0._wp … … 848 1038 849 1039 ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice) 850 z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj)1040 z_dqsb(ji,jj,jl) = rhoa(ji,jj) * rCp_air * Ch_ice(ji,jj) * wndm_ice(ji,jj) 851 1041 852 1042 ! ----------------------------! … … 863 1053 END DO 864 1054 ! 865 tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! total precipitation [kg/m2/s]866 sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! solid precipitation [kg/m2/s]867 CALL iom_put( 'snowpre', sprecip ) 868 CALL iom_put( 'precip' , tprecip ) 1055 tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1) ! total precipitation [kg/m2/s] 1056 sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1) ! solid precipitation [kg/m2/s] 1057 CALL iom_put( 'snowpre', sprecip ) ! Snow precipitation 1058 CALL iom_put( 'precip' , tprecip ) ! Total precipitation 869 1059 870 1060 ! --- evaporation --- ! … … 883 1073 ! --- heat flux associated with emp --- ! 884 1074 qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp & ! evap at sst 885 & + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp& ! liquid precip at Tair1075 & + ( tprecip(:,:) - sprecip(:,:) ) * ( ptair(:,:) - rt0 ) * rcp & ! liquid precip at Tair 886 1076 & + sprecip(:,:) * ( 1._wp - zsnw ) * & ! solid precip at min(Tair,Tsnow) 887 & ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )1077 & ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 888 1078 qemp_ice(:,:) = sprecip(:,:) * zsnw * & ! solid precip (only) 889 & ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )1079 & ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 890 1080 891 1081 ! --- total solar and non solar fluxes --- ! … … 895 1085 896 1086 ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- ! 897 qprec_ice(:,:) = rhos * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )1087 qprec_ice(:,:) = rhos * ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus ) 898 1088 899 1089 ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) --- 900 1090 DO jl = 1, jpl 901 1091 qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * rcpi * tmask(:,:,1) ) 902 ! ! But we do not have Tice => consider it at 0degC => evap=0 1092 ! ! But we do not have Tice => consider it at 0degC => evap=0 903 1093 END DO 904 1094 … … 907 1097 zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) ! zfr2 such that zfr1 + zfr2 to equal 1 908 1098 ! 909 WHERE ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) < 0.1_wp ) ! linear decrease from hi=0 to 10cm 1099 WHERE ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) < 0.1_wp ) ! linear decrease from hi=0 to 10cm 910 1100 qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) ) 911 1101 ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp ) ! constant (zfr1) when hi>10cm 912 1102 qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1 913 1103 ELSEWHERE ! zero when hs>0 914 qtr_ice_top(:,:,:) = 0._wp 1104 qtr_ice_top(:,:,:) = 0._wp 915 1105 END WHERE 916 1106 ! … … 944 1134 ENDIF 945 1135 ! 946 END SUBROUTINE blk_ice_ flx947 1136 END SUBROUTINE blk_ice_2 1137 948 1138 949 1139 SUBROUTINE blk_ice_qcn( ld_virtual_itd, ptsu, ptb, phs, phi ) … … 954 1144 !! to force sea ice / snow thermodynamics 955 1145 !! in the case conduction flux is emulated 956 !! 1146 !! 957 1147 !! ** Method : compute surface energy balance assuming neglecting heat storage 958 1148 !! following the 0-layer Semtner (1976) approach … … 979 1169 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zgfac ! enhanced conduction factor 980 1170 !!--------------------------------------------------------------------- 981 1171 982 1172 ! -------------------------------------! 983 1173 ! I Enhanced conduction factor ! … … 987 1177 ! 988 1178 zgfac(:,:,:) = 1._wp 989 1179 990 1180 IF( ld_virtual_itd ) THEN 991 1181 ! … … 993 1183 zfac2 = EXP(1._wp) * 0.5_wp * zepsilon 994 1184 zfac3 = 2._wp / zepsilon 995 ! 996 DO jl = 1, jpl 1185 ! 1186 DO jl = 1, jpl 997 1187 DO jj = 1 , jpj 998 1188 DO ji = 1, jpi … … 1002 1192 END DO 1003 1193 END DO 1004 ! 1005 ENDIF 1006 1194 ! 1195 ENDIF 1196 1007 1197 ! -------------------------------------------------------------! 1008 1198 ! II Surface temperature and conduction flux ! … … 1014 1204 DO jj = 1 , jpj 1015 1205 DO ji = 1, jpi 1016 ! 1206 ! 1017 1207 zkeff_h = zfac * zgfac(ji,jj,jl) / & ! Effective conductivity of the snow-ice system divided by thickness 1018 1208 & ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) ) … … 1031 1221 qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) 1032 1222 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) ) & 1033 1223 & * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) ) 1034 1224 1035 1225 ! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- ! 1036 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) 1226 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) 1037 1227 1038 1228 END DO 1039 1229 END DO 1040 1230 ! 1041 END DO 1042 ! 1231 END DO 1232 ! 1043 1233 END SUBROUTINE blk_ice_qcn 1044 1045 1046 SUBROUTINE Cdn10_Lupkes2012( Cd )1234 1235 1236 SUBROUTINE Cdn10_Lupkes2012( pcd ) 1047 1237 !!---------------------------------------------------------------------- 1048 1238 !! *** ROUTINE Cdn10_Lupkes2012 *** 1049 1239 !! 1050 !! ** Purpose : Recompute the neutral air-ice drag referenced at 10m 1240 !! ** Purpose : Recompute the neutral air-ice drag referenced at 10m 1051 1241 !! to make it dependent on edges at leads, melt ponds and flows. 1052 1242 !! After some approximations, this can be resumed to a dependency 1053 1243 !! on ice concentration. 1054 !! 1244 !! 1055 1245 !! ** Method : The parameterization is taken from Lupkes et al. (2012) eq.(50) 1056 1246 !! with the highest level of approximation: level4, eq.(59) … … 1064 1254 !! 1065 1255 !! This new drag has a parabolic shape (as a function of A) starting at 1066 !! Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5 1256 !! Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5 1067 1257 !! and going down to Cdi(say 1.4e-3) for A=1 1068 1258 !! … … 1074 1264 !! 1075 1265 !!---------------------------------------------------------------------- 1076 REAL(wp), DIMENSION(:,:), INTENT(inout) :: Cd1266 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pcd 1077 1267 REAL(wp), PARAMETER :: zCe = 2.23e-03_wp 1078 1268 REAL(wp), PARAMETER :: znu = 1._wp … … 1089 1279 1090 1280 ! ice-atm drag 1091 Cd(:,:) =Cd_ice + & ! pure ice drag1092 & zCe * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp) ! change due to sea-ice morphology1093 1281 pcd(:,:) = rCd_ice + & ! pure ice drag 1282 & zCe * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp) ! change due to sea-ice morphology 1283 1094 1284 END SUBROUTINE Cdn10_Lupkes2012 1095 1285 1096 1286 1097 SUBROUTINE Cdn10_Lupkes2015( Cd, Ch )1287 SUBROUTINE Cdn10_Lupkes2015( ptm_su, pslp, pcd, pch ) 1098 1288 !!---------------------------------------------------------------------- 1099 1289 !! *** ROUTINE Cdn10_Lupkes2015 *** 1100 1290 !! 1101 1291 !! ** pUrpose : Alternative turbulent transfert coefficients formulation 1102 !! between sea-ice and atmosphere with distinct momentum 1103 !! and heat coefficients depending on sea-ice concentration 1292 !! between sea-ice and atmosphere with distinct momentum 1293 !! and heat coefficients depending on sea-ice concentration 1104 1294 !! and atmospheric stability (no meltponds effect for now). 1105 !! 1295 !! 1106 1296 !! ** Method : The parameterization is adapted from Lupkes et al. (2015) 1107 1297 !! and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme, 1108 1298 !! it considers specific skin and form drags (Andreas et al. 2010) 1109 !! to compute neutral transfert coefficients for both heat and 1299 !! to compute neutral transfert coefficients for both heat and 1110 1300 !! momemtum fluxes. Atmospheric stability effect on transfert 1111 1301 !! coefficient is also taken into account following Louis (1979). … … 1116 1306 !!---------------------------------------------------------------------- 1117 1307 ! 1118 REAL(wp), DIMENSION(:,:), INTENT(inout) :: Cd 1119 REAL(wp), DIMENSION(:,:), INTENT(inout) :: Ch 1120 REAL(wp), DIMENSION(jpi,jpj) :: ztm_su, zst, zqo_sat, zqi_sat 1308 REAL(wp), DIMENSION(:,:), INTENT(in ) :: ptm_su ! sea-ice surface temperature [K] 1309 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pslp ! sea-level pressure [Pa] 1310 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pcd ! momentum transfert coefficient 1311 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pch ! heat transfert coefficient 1312 REAL(wp), DIMENSION(jpi,jpj) :: zst, zqo_sat, zqi_sat 1121 1313 ! 1122 1314 ! ECHAM6 constants … … 1146 1338 !!---------------------------------------------------------------------- 1147 1339 1148 ! mean temperature1149 WHERE( at_i_b(:,:) > 1.e-20 ) ; ztm_su(:,:) = SUM( t_su(:,:,:) * a_i_b(:,:,:) , dim=3 ) / at_i_b(:,:)1150 ELSEWHERE ; ztm_su(:,:) = rt01151 ENDWHERE1152 1153 1340 ! Momentum Neutral Transfert Coefficients (should be a constant) 1154 1341 zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2 ! Eq. 40 1155 1342 zCdn_skin_ice = ( vkarmn / LOG( rn_zu / z0_skin_ice + 1._wp ) )**2 ! Eq. 7 1156 zCdn_ice = zCdn_skin_ice ! Eq. 7 (cf Lupkes email for details)1343 zCdn_ice = zCdn_skin_ice ! Eq. 7 1157 1344 !zCdn_ice = 1.89e-3 ! old ECHAM5 value (cf Eq. 32) 1158 1345 1159 1346 ! Heat Neutral Transfert Coefficients 1160 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)1161 1347 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 1348 1162 1349 ! Atmospheric and Surface Variables 1163 1350 zst(:,:) = sst_m(:,:) + rt0 ! convert SST from Celcius to Kelvin 1164 zqo_sat(:,:) = 0.98_wp * q_sat( zst(:,:) , sf(jp_slp)%fnow(:,:,1) )! saturation humidity over ocean [kg/kg]1165 zqi_sat(:,:) = 0.98_wp * q_sat( ztm_su(:,:), sf(jp_slp)%fnow(:,:,1) )! saturation humidity over ice [kg/kg]1351 zqo_sat(:,:) = rdct_qsat_salt * q_sat( zst(:,:) , pslp(:,:) ) ! saturation humidity over ocean [kg/kg] 1352 zqi_sat(:,:) = q_sat( ptm_su(:,:), pslp(:,:) ) ! saturation humidity over ice [kg/kg] 1166 1353 ! 1167 1354 DO jj = 2, jpjm1 ! reduced loop is necessary for reproducibility … … 1169 1356 ! Virtual potential temperature [K] 1170 1357 zthetav_os = zst(ji,jj) * ( 1._wp + rctv0 * zqo_sat(ji,jj) ) ! over ocean 1171 zthetav_is = ztm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) ) ! ocean ice1358 zthetav_is = ptm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) ) ! ocean ice 1172 1359 zthetav_zu = t_zu (ji,jj) * ( 1._wp + rctv0 * q_zu(ji,jj) ) ! at zu 1173 1360 1174 1361 ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead) 1175 1362 zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj) )**2 ! over ocean 1176 1363 zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2 ! over ice 1177 1364 1178 1365 ! Momentum and Heat Neutral Transfert Coefficients 1179 1366 zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta ! Eq. 40 1180 zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) ) ! Eq. 53 1181 1182 ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead )1367 zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) ) ! Eq. 53 1368 1369 ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead ?) 1183 1370 z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water 1184 z0i = z0_skin_ice ! over ice (cf Lupkes email for details)1371 z0i = z0_skin_ice ! over ice 1185 1372 IF( zrib_o <= 0._wp ) THEN 1186 1373 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 … … 1191 1378 zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 28 1192 1379 ENDIF 1193 1380 1194 1381 IF( zrib_i <= 0._wp ) THEN 1195 1382 zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 9 … … 1199 1386 zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 27 1200 1387 ENDIF 1201 1388 1202 1389 ! Momentum Transfert Coefficients (Eq. 38) 1203 Cd(ji,jj) = zCdn_skin_ice * zfmi + &1390 pcd(ji,jj) = zCdn_skin_ice * zfmi + & 1204 1391 & 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) ) 1205 1392 1206 1393 ! Heat Transfert Coefficients (Eq. 49) 1207 Ch(ji,jj) = zChn_skin_ice * zfhi + &1394 pch(ji,jj) = zChn_skin_ice * zfhi + & 1208 1395 & 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) ) 1209 1396 ! 1210 1397 END DO 1211 1398 END DO 1212 CALL lbc_lnk_multi( 'sbcblk', Cd, 'T', 1., Ch, 'T', 1. )1399 CALL lbc_lnk_multi( 'sbcblk', pcd, 'T', 1., pch, 'T', 1. ) 1213 1400 ! 1214 1401 END SUBROUTINE Cdn10_Lupkes2015
Note: See TracChangeset
for help on using the changeset viewer.
