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