- Timestamp:
- 2019-11-29T16:59:07+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbcblk.F90
r11586 r12015 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 20 !! 4.0 ! 2019-03 (F. Lemarié & G. Samson) add ABL compatibility (ln_abl=TRUE) 21 21 !!---------------------------------------------------------------------- … … 24 24 !! sbc_blk_init : initialisation of the chosen bulk formulation as ocean surface boundary condition 25 25 !! sbc_blk : bulk formulation as ocean surface boundary condition 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 !! rho_air : density of (moist) air (depends on T_air, q_air and SLP 29 !! cp_air : specific heat of (moist) air (depends spec. hum. q_air) 30 !! q_sat : saturation humidity as a function of SLP and temperature 31 !! L_vap : latent heat of vaporization of water as a function of temperature 32 !! sea-ice case only : 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 : 33 29 !! blk_ice_1 : provide the air-ice stress 34 30 !! blk_ice_2 : provide the heat and mass fluxes at air-ice interface 35 31 !! blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux) 36 32 !! Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag 37 !! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 33 !! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 38 34 !!---------------------------------------------------------------------- 39 35 USE oce ! ocean dynamics and tracers … … 51 47 USE icethd_dh ! for CALL ice_thd_snwblow 52 48 #endif 53 USE sbcblk_algo_ncar ! => turb_ncar : NCAR - CORE (Large & Yeager, 2009) 54 USE sbcblk_algo_coare ! => turb_coare : COAREv3.0 (Fairall et al. 2003)55 USE sbcblk_algo_coare3p 5 ! => turb_coare3p5 : COAREv3.5 (Edson et al. 2013)56 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) 57 53 ! 58 54 USE iom ! I/O manager library … … 62 58 USE prtctl ! Print control 63 59 60 USE sbcblk_phy !LB: all thermodynamics functions in the marine boundary layer, rho_air, q_sat, etc... 61 62 64 63 IMPLICIT NONE 65 64 PRIVATE … … 69 68 PUBLIC blk_oce_1 ! called in sbcabl 70 69 PUBLIC blk_oce_2 ! called in sbcabl 71 PUBLIC rho_air ! called in ablmod72 PUBLIC cp_air ! called in ablmod73 70 #if defined key_si3 74 71 PUBLIC blk_ice_1 ! routine called in icesbc 75 72 PUBLIC blk_ice_2 ! routine called in icesbc 76 73 PUBLIC blk_ice_qcn ! routine called in icesbc 77 #endif 78 79 INTERFACE cp_air 80 MODULE PROCEDURE cp_air_0d, cp_air_2d 81 END INTERFACE 82 83 ! !!* Namelist namsbc_blk : bulk parameters 84 LOGICAL :: ln_NCAR ! "NCAR" algorithm (Large and Yeager 2008) 85 LOGICAL :: ln_COARE_3p0 ! "COARE 3.0" algorithm (Fairall et al. 2003) 86 LOGICAL :: ln_COARE_3p5 ! "COARE 3.5" algorithm (Edson et al. 2013) 87 LOGICAL :: ln_ECMWF ! "ECMWF" algorithm (IFS cycle 31) 88 ! ! 89 REAL(wp) :: rn_pfac ! multiplication factor for precipitation 90 REAL(wp), PUBLIC :: rn_efac !: multiplication factor for evaporation 91 REAL(wp), PUBLIC :: rn_vfac !: multiplication factor for ice/ocean velocity in the calculation of wind stress 92 REAL(wp) :: rn_zqt ! z(q,t) : height of humidity and temperature measurements 93 REAL(wp) :: rn_zu ! z(u) : height of wind measurements 94 ! ! 95 LOGICAL :: ln_Cd_L12 ! ice-atm drag = F( ice concentration ) (Lupkes et al. JGR2012) 96 LOGICAL :: ln_Cd_L15 ! ice-atm drag = F( ice concentration, atmospheric stability ) (Lupkes et al. JGR2015) 97 98 INTEGER :: nblk ! choice of the bulk algorithm 99 ! ! associated indices: 100 INTEGER, PARAMETER :: np_NCAR = 1 ! "NCAR" algorithm (Large and Yeager 2008) 101 INTEGER, PARAMETER :: np_COARE_3p0 = 2 ! "COARE 3.0" algorithm (Fairall et al. 2003) 102 INTEGER, PARAMETER :: np_COARE_3p5 = 3 ! "COARE 3.5" algorithm (Edson et al. 2013) 103 INTEGER, PARAMETER :: np_ECMWF = 4 ! "ECMWF" algorithm (IFS cycle 31) 104 105 ! !!! air parameters 106 REAL(wp) , PARAMETER :: Cp_dry = 1005.0 ! Specic heat of dry air, constant pressure [J/K/kg] 107 REAL(wp) , PARAMETER :: Cp_vap = 1860.0 ! Specic heat of water vapor, constant pressure [J/K/kg] 108 REAL(wp), PUBLIC, PARAMETER :: R_dry = 287.05_wp !: Specific gas constant for dry air [J/K/kg] 109 REAL(wp) , PARAMETER :: R_vap = 461.495_wp ! Specific gas constant for water vapor [J/K/kg] 110 REAL(wp) , PARAMETER :: reps0 = R_dry/R_vap ! ratio of gas constant for dry air and water vapor => ~ 0.622 111 REAL(wp), PUBLIC, PARAMETER :: rctv0 = R_vap/R_dry !: for virtual temperature (== (1-eps)/eps) => ~ 0.608 112 ! !!! Bulk parameters 113 REAL(wp) , PARAMETER :: cpa = 1000.5 ! specific heat of air (only used for ice fluxes now...) 114 REAL(wp) , PARAMETER :: Ls = 2.839e6 ! latent heat of sublimation 115 REAL(wp) , PARAMETER :: Stef = 5.67e-8 ! Stefan Boltzmann constant 116 REAL(wp) , PARAMETER :: rcdice = 1.4e-3 ! transfer coefficient over ice 117 REAL(wp) , PARAMETER :: albo = 0.066 ! ocean albedo assumed to be constant 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 91 ! !!* Namelist namsbc_blk : bulk parameters 92 LOGICAL :: ln_NCAR ! "NCAR" algorithm (Large and Yeager 2008) 93 LOGICAL :: ln_COARE_3p0 ! "COARE 3.0" algorithm (Fairall et al. 2003) 94 LOGICAL :: ln_COARE_3p6 ! "COARE 3.6" algorithm (Edson et al. 2013) 95 LOGICAL :: ln_ECMWF ! "ECMWF" algorithm (IFS cycle 45r1) 96 ! 97 REAL(wp) :: rn_pfac ! multiplication factor for precipitation 98 REAL(wp), PUBLIC :: rn_efac ! multiplication factor for evaporation 99 REAL(wp), PUBLIC :: rn_vfac ! multiplication factor for ice/ocean velocity in the calculation of wind stress 100 REAL(wp) :: rn_zqt ! z(q,t) : height of humidity and temperature measurements 101 REAL(wp) :: rn_zu ! z(u) : height of wind measurements 102 ! 103 LOGICAL :: ln_Cd_L12 ! ice-atm drag = F( ice concentration ) (Lupkes et al. JGR2012) 104 LOGICAL :: ln_Cd_L15 ! ice-atm drag = F( ice concentration, atmospheric stability ) (Lupkes et al. JGR2015) 118 105 ! 119 106 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: Cd_ice , Ch_ice , Ce_ice ! transfert coefficients over ice … … 121 108 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: t_zu, q_zu ! air temp. and spec. hum. at wind speed height (L15 bulk scheme) 122 109 123 !INTEGER , PUBLIC, PARAMETER :: jpfld =11 !: maximum number of files to read 124 INTEGER , PUBLIC :: jpfld !: maximum number of files to read 125 INTEGER , PUBLIC, PARAMETER :: jp_wndi = 1 !: index of 10m wind velocity (i-component) (m/s) at T-point 126 INTEGER , PUBLIC, PARAMETER :: jp_wndj = 2 !: index of 10m wind velocity (j-component) (m/s) at T-point 127 INTEGER , PUBLIC, PARAMETER :: jp_tair = 3 !: index of 10m air temperature (Kelvin) 128 INTEGER , PUBLIC, PARAMETER :: jp_humi = 4 !: index of specific humidity ( % ) 129 INTEGER , PUBLIC, PARAMETER :: jp_qsr = 5 !: index of solar heat (W/m2) 130 INTEGER , PUBLIC, PARAMETER :: jp_qlw = 6 !: index of Long wave (W/m2) 131 INTEGER , PUBLIC, PARAMETER :: jp_prec = 7 !: index of total precipitation (rain+snow) (Kg/m2/s) 132 INTEGER , PUBLIC, PARAMETER :: jp_snow = 8 !: index of snow (solid prcipitation) (kg/m2/s) 133 INTEGER , PUBLIC, PARAMETER :: jp_slp = 9 !: index of sea level pressure (Pa) 134 INTEGER , PUBLIC, PARAMETER :: jp_hpgi =10 !: index of ABL geostrophic wind or hpg (i-component) (m/s) at T-point 135 INTEGER , PUBLIC, PARAMETER :: jp_hpgj =11 !: index of ABL geostrophic wind or hpg (j-component) (m/s) at T-point 136 137 TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) :: sf ! structure of input atmospheric fields (file informations, fields read) 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 121 122 INTEGER :: nblk ! choice of the bulk algorithm 123 ! ! associated indices: 124 INTEGER, PARAMETER :: np_NCAR = 1 ! "NCAR" algorithm (Large and Yeager 2008) 125 INTEGER, PARAMETER :: np_COARE_3p0 = 2 ! "COARE 3.0" algorithm (Fairall et al. 2003) 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) 138 128 139 129 !! * Substitutions … … 165 155 !! ** Purpose : choose and initialize a bulk formulae formulation 166 156 !! 167 !! ** Method : 157 !! ** Method : 168 158 !! 169 159 !!---------------------------------------------------------------------- … … 172 162 !! 173 163 CHARACTER(len=100) :: cn_dir ! Root directory for location of atmospheric forcing files 174 !TYPE(FLD_N), DIMENSION(jpfld) :: slf_i ! array of namelist informations on the fields to read 175 TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) :: slf_i ! array of namelist informations on the fields to read 164 TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) :: slf_i ! array of namelist informations on the fields to read 176 165 TYPE(FLD_N) :: sn_wndi, sn_wndj, sn_humi, sn_qsr ! informations about the fields to be read 177 166 TYPE(FLD_N) :: sn_qlw , sn_tair, sn_prec, sn_snow ! " " … … 179 168 NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw , & ! input fields 180 169 & sn_tair, sn_prec, sn_snow, sn_slp, sn_hpgi, sn_hpgj, & 181 & ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF, & ! bulk algorithm 182 & cn_dir , rn_zqt, rn_zu, & 183 & rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15 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 184 174 !!--------------------------------------------------------------------- 185 175 ! … … 187 177 IF( sbc_blk_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' ) 188 178 ! 189 ! !** read bulk namelist 179 ! !** read bulk namelist 190 180 REWIND( numnam_ref ) !* Namelist namsbc_blk in reference namelist : bulk parameters 191 181 READ ( numnam_ref, namsbc_blk, IOSTAT = ios, ERR = 901) … … 200 190 ! !** initialization of the chosen bulk formulae (+ check) 201 191 ! !* select the bulk chosen in the namelist and check the choice 202 ioptio = 0 203 IF( ln_NCAR ) THEN ; nblk = np_NCAR ; ioptio = ioptio + 1 ; ENDIF 204 IF( ln_COARE_3p0 ) THEN ; nblk = np_COARE_3p0 ; ioptio = ioptio + 1 ; ENDIF 205 IF( ln_COARE_3p5 ) THEN ; nblk = np_COARE_3p5 ; ioptio = ioptio + 1 ; ENDIF 206 IF( ln_ECMWF ) THEN ; nblk = np_ECMWF ; ioptio = ioptio + 1 ; ENDIF 207 ! 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 208 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_NCAR .AND. (ln_skin_cs .OR. ln_skin_wl) ) & 209 & CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with NCAR algorithm!' ) 210 ! 211 ioptio = 0 212 IF( ln_humi_sph ) THEN 213 nhumi = np_humi_sph ; ioptio = ioptio + 1 214 ENDIF 215 IF( ln_humi_dpt ) THEN 216 nhumi = np_humi_dpt ; ioptio = ioptio + 1 217 ENDIF 218 IF( ln_humi_rlh ) THEN 219 nhumi = np_humi_rlh ; ioptio = ioptio + 1 220 ENDIF 221 IF( ioptio /= 1 ) CALL ctl_stop( 'sbc_blk_init: Choose one and only one type of air humidity' ) 209 222 ! 210 223 IF( ln_dm2dc ) THEN !* check: diurnal cycle on Qsr 211 224 IF( sn_qsr%freqh /= 24. ) CALL ctl_stop( 'sbc_blk_init: ln_dm2dc=T only with daily short-wave input' ) 212 IF( sn_qsr%ln_tint ) THEN 225 IF( sn_qsr%ln_tint ) THEN 213 226 CALL ctl_warn( 'sbc_blk_init: ln_dm2dc=T daily qsr time interpolation done by sbcdcy module', & 214 227 & ' ==> We force time interpolation = .false. for qsr' ) … … 234 247 ALLOCATE( sf(jpfld), STAT=ierror ) 235 248 IF( ierror > 0 ) CALL ctl_stop( 'STOP', 'sbc_blk_init: unable to allocate sf structure' ) 236 ! !- fill the bulk structure with namelist informations237 CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' )238 249 ! 239 250 DO jfpr= 1, jpfld … … 258 269 ENDIF 259 270 END DO 260 ! 261 IF( ln_wave ) THEN ! surface waves 262 IF( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) ) & ! Activated wave module but neither drag nor stokes drift activated 263 & CALL ctl_stop( 'sbc_blk_init: Ask for wave coupling but ln_cdgw=ln_sdw=ln_tauwoc=ln_stcor=F' ) 264 IF( ln_cdgw .AND. .NOT.ln_NCAR ) & ! drag coefficient read from wave model only with NCAR bulk formulae 265 & CALL ctl_stop( 'sbc_blk_init: drag coefficient read from wave model need NCAR bulk formulae') 266 IF( ln_stcor .AND. .NOT.ln_sdw ) & 267 CALL ctl_stop( 'sbc_blk_init: Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T') 271 ! !- fill the bulk structure with namelist informations 272 CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' ) 273 ! 274 IF( ln_wave ) THEN 275 !Activated wave module but neither drag nor stokes drift activated 276 IF( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) ) THEN 277 CALL ctl_stop( 'STOP', 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauwoc=F, ln_stcor=F' ) 278 !drag coefficient read from wave model definable only with mfs bulk formulae and core 279 ELSEIF(ln_cdgw .AND. .NOT. ln_NCAR ) THEN 280 CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR and CORE bulk formulae') 281 ELSEIF(ln_stcor .AND. .NOT. ln_sdw) THEN 282 CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T') 283 ENDIF 268 284 ELSE 269 IF( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) THEN 270 CALL ctl_warn( 'sbc_blk_init: ln_wave=F, set all wave-related namelist parameter to FALSE') 271 ln_cdgw =.FALSE. ; ln_sdw =.FALSE. ; ln_tauwoc =.FALSE. ; ln_stcor =.FALSE. 272 ENDIF 273 ENDIF 285 IF( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) & 286 & CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ', & 287 & 'with drag coefficient (ln_cdgw =T) ' , & 288 & 'or Stokes Drift (ln_sdw=T) ' , & 289 & 'or ocean stress modification due to waves (ln_tauwoc=T) ', & 290 & 'or Stokes-Coriolis term (ln_stcori=T)' ) 291 ENDIF 274 292 ! 275 293 IF( ln_abl ) THEN ! ABL: read 3D fields for wind, temperature, humidity and pressure gradient … … 281 299 ! 282 300 ! set transfer coefficients to default sea-ice values 283 Cd_ice(:,:) = r cdice284 Ch_ice(:,:) = r cdice285 Ce_ice(:,:) = r cdice301 Cd_ice(:,:) = rCd_ice 302 Ch_ice(:,:) = rCd_ice 303 Ce_ice(:,:) = rCd_ice 286 304 ! 287 305 IF(lwp) THEN !** Control print 288 306 ! 289 WRITE(numout,*) !* namelist 307 WRITE(numout,*) !* namelist 290 308 WRITE(numout,*) ' Namelist namsbc_blk (other than data information):' 291 309 WRITE(numout,*) ' "NCAR" algorithm (Large and Yeager 2008) ln_NCAR = ', ln_NCAR 292 310 WRITE(numout,*) ' "COARE 3.0" algorithm (Fairall et al. 2003) ln_COARE_3p0 = ', ln_COARE_3p0 293 WRITE(numout,*) ' "COARE 3. 5" algorithm (Edson et al. 2013) ln_COARE_3p5 = ', ln_COARE_3p0294 WRITE(numout,*) ' "ECMWF" algorithm (IFS cycle 31)ln_ECMWF = ', ln_ECMWF311 WRITE(numout,*) ' "COARE 3.6" algorithm (Fairall 2018 + Edson al 2013)ln_COARE_3p6 = ', ln_COARE_3p6 312 WRITE(numout,*) ' "ECMWF" algorithm (IFS cycle 45r1) ln_ECMWF = ', ln_ECMWF 295 313 WRITE(numout,*) ' Air temperature and humidity reference height (m) rn_zqt = ', rn_zqt 296 314 WRITE(numout,*) ' Wind vector reference height (m) rn_zu = ', rn_zu … … 306 324 CASE( np_NCAR ) ; WRITE(numout,*) ' ==>>> "NCAR" algorithm (Large and Yeager 2008)' 307 325 CASE( np_COARE_3p0 ) ; WRITE(numout,*) ' ==>>> "COARE 3.0" algorithm (Fairall et al. 2003)' 308 CASE( np_COARE_3p5 ) ; WRITE(numout,*) ' ==>>> "COARE 3.5" algorithm (Edson et al. 2013)' 309 CASE( np_ECMWF ) ; WRITE(numout,*) ' ==>>> "ECMWF" algorithm (IFS cycle 31)' 326 CASE( np_COARE_3p6 ) ; WRITE(numout,*) ' ==>>> "COARE 3.6" algorithm (Fairall 2018+Edson et al. 2013)' 327 CASE( np_ECMWF ) ; WRITE(numout,*) ' ==>>> "ECMWF" algorithm (IFS cycle 45r1)' 328 END SELECT 329 ! 330 WRITE(numout,*) 331 WRITE(numout,*) ' use cool-skin parameterization (SSST) ln_skin_cs = ', ln_skin_cs 332 WRITE(numout,*) ' use warm-layer parameterization (SSST) ln_skin_wl = ', ln_skin_wl 333 ! 334 WRITE(numout,*) 335 SELECT CASE( nhumi ) !* Print the choice of air humidity 336 CASE( np_humi_sph ) ; WRITE(numout,*) ' ==>>> air humidity is SPECIFIC HUMIDITY [kg/kg]' 337 CASE( np_humi_dpt ) ; WRITE(numout,*) ' ==>>> air humidity is DEW-POINT TEMPERATURE [K]' 338 CASE( np_humi_rlh ) ; WRITE(numout,*) ' ==>>> air humidity is RELATIVE HUMIDITY [%]' 310 339 END SELECT 311 340 ! … … 355 384 CALL fld_read( kt, nn_fsbc, sf ) ! input fields provided at the current time-step 356 385 ! 357 ! ! compute the surface ocean fluxes using bulk formulea358 IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN359 !386 IF( kt == nit000 ) tsk(:,:) = sst_m(:,:)*tmask(:,:,1) ! no previous estimate of skin temperature => using bulk SST 387 ! 388 IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 360 389 CALL blk_oce_1( kt, sf(jp_wndi)%fnow(:,:,1), sf(jp_wndj)%fnow(:,:,1), & ! <<= in 361 390 & sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), & ! <<= in 362 391 & sf(jp_slp )%fnow(:,:,1), sst_m, ssu_m, ssv_m, & ! <<= in 392 & sf(jp_qsr )%fnow(:,:,1), sf(jp_qlw )%fnow(:,:,1), & ! <<= in (wl/cs) 363 393 & zssq, zcd_du, zsen, zevp ) ! =>> out 364 394 … … 368 398 & zsen, zevp ) ! <=> in out 369 399 ENDIF 370 400 ! 371 401 #if defined key_cice 372 IF( MOD( kt -1, nn_fsbc ) == 0 ) THEN402 IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 373 403 qlw_ice(:,:,1) = sf(jp_qlw )%fnow(:,:,1) 374 IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) 375 ELSE ; qsr_ice(:,:,1) = sf(jp_qsr)%fnow(:,:,1) 376 ENDIF 404 IF( ln_dm2dc ) THEN 405 qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) 406 ELSE 407 qsr_ice(:,:,1) = sf(jp_qsr)%fnow(:,:,1) 408 ENDIF 377 409 tatm_ice(:,:) = sf(jp_tair)%fnow(:,:,1) 378 qatm_ice(:,:) = sf(jp_humi)%fnow(:,:,1) 410 411 SELECT CASE( nhumi ) 412 CASE( np_humi_sph ) 413 qatm_ice(:,:) = sf(jp_humi)%fnow(:,:,1) 414 CASE( np_humi_dpt ) 415 qatm_ice(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 416 CASE( np_humi_rlh ) 417 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 418 END SELECT 419 379 420 tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac 380 421 sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac … … 387 428 388 429 389 SUBROUTINE blk_oce_1( kt, pwndi, pwndj, ptair, phumi, & ! inp 390 & pslp, pst , pu , pv, & ! inp 391 & pssq, pcd_du, psen, pevp ) ! out 430 SUBROUTINE blk_oce_1( kt, pwndi, pwndj , ptair, phumi, & ! inp 431 & pslp , pst , pu , pv, & ! inp 432 & pqsr , pqlw , & ! inp 433 & pssq , pcd_du, psen , pevp ) ! out 392 434 !!--------------------------------------------------------------------- 393 435 !! *** ROUTINE blk_oce_1 *** 394 436 !! 395 !! ** Purpose : Computes Cd x |U|, Ch x |U|, Ce x |U| for ABL integration 396 !! 397 !! ** Method : bulk formulae using atmospheric 398 !! fields from the ABL model at previous time-step 437 !! ** Purpose : if ln_blk=T, computes surface momentum, heat and freshwater fluxes 438 !! if ln_abl=T, computes Cd x |U|, Ch x |U|, Ce x |U| for ABL integration 439 !! 440 !! ** Method : bulk formulae using atmospheric fields from : 441 !! if ln_blk=T, atmospheric fields read in sbc_read 442 !! if ln_abl=T, the ABL model at previous time-step 399 443 !! 400 444 !! ** Outputs : - pssq : surface humidity used to compute latent heat flux (kg/kg) … … 412 456 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pu ! surface current at U-point (i-component) [m/s] 413 457 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pv ! surface current at V-point (j-component) [m/s] 458 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pqsr ! 459 REAL(wp), INTENT(in ), DIMENSION(:,:) :: pqlw ! 414 460 REAL(wp), INTENT( out), DIMENSION(:,:) :: pssq ! specific humidity at pst [kg/kg] 415 461 REAL(wp), INTENT( out), DIMENSION(:,:) :: pcd_du ! Cd x |dU| at T-points [m/s] … … 423 469 REAL(wp), DIMENSION(jpi,jpj) :: zU_zu ! bulk wind speed at height zu [m/s] 424 470 REAL(wp), DIMENSION(jpi,jpj) :: ztpot ! potential temperature of air at z=rn_zqt [K] 425 REAL(wp), DIMENSION(jpi,jpj) :: z rhoa ! density of air [kg/m^3]471 REAL(wp), DIMENSION(jpi,jpj) :: zqair ! specific humidity of air at z=rn_zqt [kg/kg] 426 472 REAL(wp), DIMENSION(jpi,jpj) :: zcd_oce ! momentum transfert coefficient over ocean 427 473 REAL(wp), DIMENSION(jpi,jpj) :: zch_oce ! sensible heat transfert coefficient over ocean 428 474 REAL(wp), DIMENSION(jpi,jpj) :: zce_oce ! latent heat transfert coefficient over ocean 475 REAL(wp), DIMENSION(jpi,jpj) :: zqla ! latent heat flux 476 REAL(wp), DIMENSION(jpi,jpj) :: zztmp1, zztmp2 429 477 !!--------------------------------------------------------------------- 430 478 ! … … 460 508 461 509 ! ----------------------------------------------------------------------------- ! 462 ! I Turbulent FLUXES!510 ! I Solar FLUX ! 463 511 ! ----------------------------------------------------------------------------- ! 464 512 465 ! ... specific humidity at SST and IST tmask( 466 pssq(:,:) = 0.98 * q_sat( zst(:,:), pslp(:,:) ) 467 ! 513 ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle ! Short Wave 514 zztmp = 1. - albo 515 IF( ln_dm2dc ) THEN 516 qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1) 517 ELSE 518 qsr(:,:) = zztmp * sf(jp_qsr)%fnow(:,:,1) * tmask(:,:,1) 519 ENDIF 520 521 522 ! ----------------------------------------------------------------------------- ! 523 ! II Turbulent FLUXES ! 524 ! ----------------------------------------------------------------------------- ! 525 526 ! specific humidity at SST 527 pssq(:,:) = rdct_qsat_salt * q_sat( zst(:,:), pslp(:,:) ) 528 529 IF( ln_skin_cs .OR. ln_skin_wl ) THEN 530 zztmp1(:,:) = zst(:,:) 531 zztmp2(:,:) = pssq(:,:) 532 ENDIF 533 534 ! specific humidity of air at "rn_zqt" m above the sea 535 SELECT CASE( nhumi ) 536 CASE( np_humi_sph ) 537 zqair(:,:) = phumi(:,:) ! what we read in file is already a spec. humidity! 538 CASE( np_humi_dpt ) 539 !IF(lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of d_air and slp !' !LBrm 540 zqair(:,:) = q_sat( phumi(:,:), pslp(:,:) ) 541 CASE( np_humi_rlh ) 542 !IF(lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of RH, t_air and slp !' !LBrm 543 zqair(:,:) = q_air_rh( 0.01_wp*phumi(:,:), ptair(:,:), pslp(:,:) ) !LB: 0.01 => RH is % percent in file 544 END SELECT 545 ! 546 ! potential temperature of air at "rn_zqt" m above the sea 468 547 IF( ln_abl ) THEN 469 548 ztpot = ptair(:,:) … … 472 551 ! (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2 473 552 ! (since reanalysis products provide T at z, not theta !) 474 ztpot = ptair(:,:) + gamma_moist( ptair(:,:), phumi(:,:) ) * rn_zqt 475 ENDIF 476 477 SELECT CASE( nblk ) !== transfer coefficients ==! Cd, Ch, Ce at T-point 478 ! 479 CASE( np_NCAR ) ; CALL turb_ncar ( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm, & ! NCAR-COREv2 480 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, Cdn_oce, Chn_oce, Cen_oce ) 481 CASE( np_COARE_3p0 ) ; CALL turb_coare ( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm, & ! COARE v3.0 482 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, Cdn_oce, Chn_oce, Cen_oce ) 483 CASE( np_COARE_3p5 ) ; CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm, & ! COARE v3.5 484 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, Cdn_oce, Chn_oce, Cen_oce ) 485 CASE( np_ECMWF ) ; CALL turb_ecmwf ( rn_zqt, rn_zu, zst, ztpot, pssq, phumi, wndm, & ! ECMWF 486 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, Cdn_oce, Chn_oce, Cen_oce ) 487 CASE DEFAULT 488 CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 489 END SELECT 490 491 492 IF ( ln_abl ) THEN !== ABL formulation ==! multiplication by rho_air and turbulent fluxes computation done in ablstp 493 !! FL do we need this multiplication by tmask ... ??? 494 DO jj = 1, jpj 495 DO ji = 1, jpi 496 zztmp = zU_zu(ji,jj) !* tmask(ji,jj,1) 497 wndm(ji,jj) = zztmp ! Store zU_zu in wndm to compute ustar2 in ablmod 498 pcd_du(ji,jj) = zztmp * zcd_oce(ji,jj) 499 psen(ji,jj) = zztmp * zch_oce(ji,jj) 500 pevp(ji,jj) = zztmp * zce_oce(ji,jj) 501 END DO 502 END DO 503 504 ELSE !== BLK formulation ==! turbulent fluxes computation 553 ztpot = ptair(:,:) + gamma_moist( ptair(:,:), zqair(:,:) ) * rn_zqt 554 ENDIF 555 556 557 558 IF( ln_skin_cs .OR. ln_skin_wl ) THEN 559 560 SELECT CASE( nblk ) !== transfer coefficients ==! Cd, Ch, Ce at T-point 561 562 CASE( np_COARE_3p0 ) 563 CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 564 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 565 & Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 566 567 CASE( np_COARE_3p6 ) 568 CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 569 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 570 & Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 571 572 CASE( np_ECMWF ) 573 CALL turb_ecmwf ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 574 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce, & 575 & Qsw=qsr(:,:), rad_lw=pqlw(:,:), slp=pslp(:,:) ) 576 577 CASE DEFAULT 578 CALL ctl_stop( 'STOP', 'sbc_oce: unsuported bulk formula selection for "ln_skin_*==.true."' ) 579 END SELECT 580 581 !! In the presence of sea-ice we forget about the cool-skin/warm-layer update of zst and pssq: 582 WHERE ( fr_i < 0.001_wp ) 583 ! zst and pssq have been updated by cool-skin/warm-layer scheme and we keep it!!! 584 zst(:,:) = zst(:,:)*tmask(:,:,1) 585 pssq(:,:) = pssq(:,:)*tmask(:,:,1) 586 ELSEWHERE 587 ! we forget about the update... 588 zst(:,:) = zztmp1(:,:) !LB: using what we backed up before skin-algo 589 pssq(:,:) = zztmp2(:,:) !LB: " " " 590 END WHERE 591 592 !LB: Update of tsk, the "official" array for skin temperature 593 tsk(:,:) = zst(:,:) 594 595 ELSE !IF( ln_skin_cs .OR. ln_skin_wl ) 596 597 SELECT CASE( nblk ) !== transfer coefficients ==! Cd, Ch, Ce at T-point 598 ! 599 CASE( np_NCAR ) 600 CALL turb_ncar ( rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, & 601 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 602 603 CASE( np_COARE_3p0 ) 604 CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 605 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 606 607 CASE( np_COARE_3p6 ) 608 CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 609 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 610 611 CASE( np_ECMWF ) 612 CALL turb_ecmwf ( kt, rn_zqt, rn_zu, zst, ztpot, pssq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 613 & zcd_oce, zch_oce, zce_oce, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 614 615 CASE DEFAULT 616 CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 617 END SELECT 618 619 ENDIF ! IF( ln_skin_cs .OR. ln_skin_wl ) 620 621 !! CALL iom_put( "Cd_oce", zcd_oce) ! output value of pure ocean-atm. transfer coef. 622 !! CALL iom_put( "Ch_oce", zch_oce) ! output value of pure ocean-atm. transfer coef. 623 624 IF( ABS(rn_zu - rn_zqt) < 0.1_wp ) THEN 625 !! If zu == zt, then ensuring once for all that: 626 t_zu(:,:) = ztpot(:,:) 627 q_zu(:,:) = zqair(:,:) 628 ENDIF 629 630 631 ! Turbulent fluxes over ocean => BULK_FORMULA @ sbcblk_phy.F90 632 ! ------------------------------------------------------------- 633 634 IF (ln_abl) THEN 635 636 CALL BULK_FORMULA( rn_zu, zst(:,:), pssq(:,:), t_zu(:,:), q_zu(:,:), & 637 & zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:), & 638 & wndm(:,:), zU_zu(:,:), pslp(:,:), & 639 & pcd_du(:,:), psen(:,:), pevp(:,:), & 640 & prhoa=rhoa(:,:) ) 641 642 ELSE 643 644 CALL BULK_FORMULA( rn_zu, zst(:,:), pssq(:,:), t_zu(:,:), q_zu(:,:), & 645 & zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:), & 646 & wndm(:,:), zU_zu(:,:), pslp(:,:), & 647 & taum(:,:), psen(:,:), zqla(:,:), & 648 & pEvap=pevp(:,:), prhoa=rhoa(:,:) ) 649 650 zqla(:,:) = zqla(:,:) * tmask(:,:,1) 651 psen(:,:) = psen(:,:) * tmask(:,:,1) 652 taum(:,:) = taum(:,:) * tmask(:,:,1) 653 pevp(:,:) = pevp(:,:) * tmask(:,:,1) 505 654 506 ! ! Compute true air density : 507 IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN ! At zu: (probably useless to remove zrho*grav*rn_zu from SLP...) 508 zrhoa(:,:) = rho_air( t_zu(:,:) , q_zu(:,:) , pslp(:,:) ) 509 ELSE ! At zt: 510 zrhoa(:,:) = rho_air( ptair(:,:), phumi(:,:), pslp(:,:) ) 511 END IF 512 513 DO jj = 1, jpj 514 DO ji = 1, jpi 515 zztmp = zrhoa(ji,jj) * zU_zu(ji,jj) 516 !!gm to be done by someone: check the efficiency of the call of cp_air within do loops 517 psen (ji,jj) = cp_air( q_zu(ji,jj) ) * zztmp * zch_oce(ji,jj) * ( zst(ji,jj) - t_zu(ji,jj) ) 518 pevp (ji,jj) = rn_efac*MAX( 0._wp, zztmp * zce_oce(ji,jj) * ( pssq(ji,jj) - q_zu(ji,jj) ) ) 519 zztmp = zztmp * zcd_oce(ji,jj) 520 pcd_du(ji,jj) = zztmp 521 taum (ji,jj) = zztmp * wndm (ji,jj) 522 zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 523 zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj) 524 END DO 525 END DO 655 ! Tau i and j component on T-grid points, using array "zcd_oce" as a temporary array... 656 zcd_oce = 0._wp 657 WHERE ( wndm > 0._wp ) zcd_oce = taum / wndm 658 zwnd_i = zcd_oce * zwnd_i 659 zwnd_j = zcd_oce * zwnd_j 526 660 527 661 CALL iom_put( "taum_oce", taum ) ! output wind stress module 528 662 529 663 ! ... utau, vtau at U- and V_points, resp. 530 664 ! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines … … 539 673 END DO 540 674 CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. ) 541 542 675 IF(ln_ctl) THEN 543 676 CALL prt_ctl( tab2d_1=wndm , clinfo1=' blk_oce_1: wndm : ') … … 554 687 ENDIF 555 688 ! 556 END SUBROUTINE blk_oce_1689 END SUBROUTINE blk_oce_1 557 690 558 691 … … 593 726 zst(:,:) = pst(:,:) + rt0 ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K) 594 727 728 595 729 ! ----------------------------------------------------------------------------- ! 596 ! I Radiative FLUXES!730 ! III Net longwave radiative FLUX ! 597 731 ! ----------------------------------------------------------------------------- ! 598 732 599 ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle ! Short Wave 600 zztmp = 1._wp - albo 601 IF( ln_dm2dc ) THEN ; qsr(:,:) = zztmp * sbc_dcy( pqsr(:,:) ) * tmask(:,:,1) 602 ELSE ; qsr(:,:) = zztmp * pqsr(:,:) * tmask(:,:,1) 603 ENDIF 604 605 zqlw(:,:) = ( pqlw(:,:) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:) ) * tmask(:,:,1) ! Long Wave 733 !! LB: now moved after Turbulent fluxes because must use the skin temperature rather that the SST 734 !! (zst is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.) 735 zqlw(:,:) = emiss_w * ( pqlw(:,:) - stefan*zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:) ) * tmask(:,:,1) ! Net radiative longwave flux 606 736 607 737 ! Turbulent fluxes over ocean 608 738 ! ----------------------------- 609 739 610 zqla(:,:) = L_vap( zst(:,:) ) * pevp(:,:) ! Latent Heat flux 740 zqla(:,:) = L_vap( zst(:,:) ) * pevp(:,:) ! Latent Heat flux !!GS: possibility to add a global qla to avoid recomputation after abl update 611 741 612 742 IF(ln_ctl) THEN 613 743 CALL prt_ctl( tab2d_1=zqla , clinfo1=' blk_oce_2: zqla : ' ) 614 744 CALL prt_ctl( tab2d_1=zqlw , clinfo1=' blk_oce_2: zqlw : ', tab2d_2=qsr, clinfo2=' qsr : ' ) 745 615 746 ENDIF 616 747 617 748 ! ----------------------------------------------------------------------------- ! 618 ! I IITotal FLUXES !749 ! IV Total FLUXES ! 619 750 ! ----------------------------------------------------------------------------- ! 620 751 ! 621 emp (:,:) = ( pevp(:,:) - pprec(:,:) * rn_pfac ) * tmask(:,:,1) ! mass flux (evap. - precip.) 622 ! 623 zz1 = rn_pfac * rLfus 624 zz2 = rn_pfac * rcp 625 zz3 = rn_pfac * rcpi 626 zztmp = 1._wp - albo 627 qns(:,:) = ( zqlw(:,:) - psen(:,:) - zqla(:,:) & ! Downward Non Solar 628 & - psnow(:,:) * zztmp & ! remove latent melting heat for solid precip 629 & - pevp(:,:) * pst(:,:) * rcp & ! remove evap heat content at SST 630 & + ( pprec(:,:) - psnow(:,:) ) * ( ptair(:,:) - rt0 ) * zz2 & ! add liquid precip heat content at Tair 631 & + psnow(:,:) * ( MIN( ptair(:,:), rt0 ) - rt0 ) * zz3 & ! add solid precip heat content at min(Tair,Tsnow) 632 & ) * tmask(:,:,1) 752 emp (:,:) = ( pevp(:,:) & ! mass flux (evap. - precip.) 753 & - pprec(:,:) * rn_pfac ) * tmask(:,:,1) 754 ! 755 qns(:,:) = zqlw(:,:) + psen(:,:) + zqla(:,:) & ! Downward Non Solar 756 & - psnow(:,:) * rn_pfac * rLfus & ! remove latent melting heat for solid precip 757 & - pevp(:,:) * pst(:,:) * rcp & ! remove evap heat content at SST !LB??? pst is Celsius !? 758 & + ( pprec(:,:) - psnow(:,:) ) * rn_pfac & ! add liquid precip heat content at Tair 759 & * ( ptair(:,:) - rt0 ) * rcp & 760 & + psnow(:,:) * rn_pfac & ! add solid precip heat content at min(Tair,Tsnow) 761 & * ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi 762 qns(:,:) = qns(:,:) * tmask(:,:,1) 633 763 ! 634 764 #if defined key_si3 635 qns_oce(:,:) = zqlw(:,:) - psen(:,:) -zqla(:,:) ! non solar without emp (only needed by SI3)765 qns_oce(:,:) = zqlw(:,:) + psen(:,:) + zqla(:,:) ! non solar without emp (only needed by SI3) 636 766 qsr_oce(:,:) = qsr(:,:) 637 767 #endif 638 768 ! 639 IF ( nn_ice == 0 ) THEN 640 CALL iom_put( "qlw_oce" , zqlw ) ! output downward longwave heat over the ocean 641 CALL iom_put( "qsb_oce" , - psen ) ! output downward sensible heat over the ocean 642 CALL iom_put( "qla_oce" , - zqla ) ! output downward latent heat over the ocean 643 CALL iom_put( "qemp_oce", qns-zqlw+psen+zqla ) ! output downward heat content of E-P over the ocean 644 CALL iom_put( "qns_oce" , qns ) ! output downward non solar heat over the ocean 645 CALL iom_put( "qsr_oce" , qsr ) ! output downward solar heat over the ocean 646 CALL iom_put( "qt_oce" , qns+qsr ) ! output total downward heat over the ocean 647 tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s] 648 sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s] 649 CALL iom_put( 'snowpre', sprecip ) ! Snow 650 CALL iom_put( 'precip' , tprecip ) ! Total precipitation 769 CALL iom_put( "rho_air" , rhoa ) ! output air density (kg/m^3) !#LB 770 CALL iom_put( "qlw_oce" , zqlw ) ! output downward longwave heat over the ocean 771 CALL iom_put( "qsb_oce" , psen ) ! output downward sensible heat over the ocean 772 CALL iom_put( "qla_oce" , zqla ) ! output downward latent heat over the ocean 773 CALL iom_put( "evap_oce" , pevp ) ! evaporation 774 CALL iom_put( "qemp_oce" , qns-zqlw-psen-zqla ) ! output downward heat content of E-P over the ocean 775 CALL iom_put( "qns_oce" , qns ) ! output downward non solar heat over the ocean 776 CALL iom_put( "qsr_oce" , qsr ) ! output downward solar heat over the ocean 777 CALL iom_put( "qt_oce" , qns+qsr ) ! output total downward heat over the ocean 778 tprecip(:,:) = pprec(:,:) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s] 779 sprecip(:,:) = psnow(:,:) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s] 780 CALL iom_put( 'snowpre', sprecip ) ! Snow 781 CALL iom_put( 'precip' , tprecip ) ! Total precipitation 782 IF( ln_skin_cs .OR. ln_skin_wl ) THEN 783 CALL iom_put( "t_skin" , (zst - rt0) * tmask(:,:,1) ) ! T_skin in Celsius 784 CALL iom_put( "dt_skin" , (zst - pst - rt0) * tmask(:,:,1) ) ! T_skin - SST temperature difference... 651 785 ENDIF 652 786 ! … … 660 794 661 795 662 FUNCTION rho_air( ptak, pqa, pslp )663 !!-------------------------------------------------------------------------------664 !! *** FUNCTION rho_air ***665 !!666 !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere667 !!668 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)669 !!-------------------------------------------------------------------------------670 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K]671 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg]672 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp ! pressure in [Pa]673 REAL(wp), DIMENSION(jpi,jpj) :: rho_air ! density of moist air [kg/m^3]674 !!-------------------------------------------------------------------------------675 !676 rho_air = pslp / ( R_dry*ptak * ( 1._wp + rctv0*pqa ) )677 !678 END FUNCTION rho_air679 680 681 FUNCTION cp_air_0d( pqa )682 !!-------------------------------------------------------------------------------683 !! *** FUNCTION cp_air ***684 !!685 !! ** Purpose : Compute specific heat (Cp) of moist air686 !!687 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)688 !!-------------------------------------------------------------------------------689 REAL(wp), INTENT(in) :: pqa ! air specific humidity [kg/kg]690 REAL(wp) :: cp_air_0d! specific heat of moist air [J/K/kg]691 !!-------------------------------------------------------------------------------692 !693 Cp_air_0d = Cp_dry + Cp_vap * pqa694 !695 END FUNCTION cp_air_0d696 697 698 FUNCTION cp_air_2d( pqa )699 !!-------------------------------------------------------------------------------700 !! *** FUNCTION cp_air ***701 !!702 !! ** Purpose : Compute specific heat (Cp) of moist air703 !!704 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)705 !!-------------------------------------------------------------------------------706 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg]707 REAL(wp), DIMENSION(jpi,jpj) :: cp_air_2d! specific heat of moist air [J/K/kg]708 !!-------------------------------------------------------------------------------709 !710 Cp_air_2d = Cp_dry + Cp_vap * pqa711 !712 END FUNCTION cp_air_2d713 714 715 FUNCTION q_sat( ptak, pslp )716 !!----------------------------------------------------------------------------------717 !! *** FUNCTION q_sat ***718 !!719 !! ** Purpose : Specific humidity at saturation in [kg/kg]720 !! Based on accurate estimate of "e_sat"721 !! aka saturation water vapor (Goff, 1957)722 !!723 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)724 !!----------------------------------------------------------------------------------725 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K]726 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pslp ! sea level atmospheric pressure [Pa]727 REAL(wp), DIMENSION(jpi,jpj) :: q_sat ! Specific humidity at saturation [kg/kg]728 !729 INTEGER :: ji, jj ! dummy loop indices730 REAL(wp) :: ze_sat, ztmp ! local scalar731 !!----------------------------------------------------------------------------------732 !733 DO jj = 1, jpj734 DO ji = 1, jpi735 !736 ztmp = rt0 / ptak(ji,jj)737 !738 ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957)739 ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0) &740 & + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) ) &741 & + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614 )742 !743 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]744 !745 END DO746 END DO747 !748 END FUNCTION q_sat749 750 751 FUNCTION gamma_moist( ptak, pqa )752 !!----------------------------------------------------------------------------------753 !! *** FUNCTION gamma_moist ***754 !!755 !! ** Purpose : Compute the moist adiabatic lapse-rate.756 !! => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate757 !! => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html758 !!759 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)760 !!----------------------------------------------------------------------------------761 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K]762 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! specific humidity [kg/kg]763 REAL(wp), DIMENSION(jpi,jpj) :: gamma_moist ! moist adiabatic lapse-rate764 !765 INTEGER :: ji, jj ! dummy loop indices766 REAL(wp) :: zrv, ziRT ! local scalar767 !!----------------------------------------------------------------------------------768 !769 DO jj = 1, jpj770 DO ji = 1, jpi771 zrv = pqa(ji,jj) / (1. - pqa(ji,jj))772 ziRT = 1. / (R_dry*ptak(ji,jj)) ! 1/RT773 gamma_moist(ji,jj) = grav * ( 1. + rLevap*zrv*ziRT ) / ( Cp_dry + rLevap*rLevap*zrv*reps0*ziRT/ptak(ji,jj) )774 END DO775 END DO776 !777 END FUNCTION gamma_moist778 779 780 FUNCTION L_vap( psst )781 !!---------------------------------------------------------------------------------782 !! *** FUNCTION L_vap ***783 !!784 !! ** Purpose : Compute the latent heat of vaporization of water from temperature785 !!786 !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)787 !!----------------------------------------------------------------------------------788 REAL(wp), DIMENSION(jpi,jpj) :: L_vap ! latent heat of vaporization [J/kg]789 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst ! water temperature [K]790 !!----------------------------------------------------------------------------------791 !792 L_vap = ( 2.501 - 0.00237 * ( psst(:,:) - rt0) ) * 1.e6793 !794 END FUNCTION L_vap795 796 796 #if defined key_si3 797 797 !!---------------------------------------------------------------------- … … 802 802 !! blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux) 803 803 !! Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag 804 !! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 804 !! Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag 805 805 !!---------------------------------------------------------------------- 806 806 … … 835 835 REAL(wp) :: zootm_su ! sea-ice surface mean temperature 836 836 REAL(wp) :: zztmp1, zztmp2 ! temporary arrays 837 REAL(wp), DIMENSION(jpi,jpj) :: zrhoa ! transfer coefficient for momentum (tau)838 837 REAL(wp), DIMENSION(jpi,jpj) :: zcd_dui ! transfer coefficient for momentum (tau) 839 838 !!--------------------------------------------------------------------- 839 ! 840 840 841 841 ! ------------------------------------------------------------ ! … … 858 858 Ce_ice(:,:) = Cd_ice(:,:) 859 859 ELSEIF( ln_Cd_L15 ) THEN ! calculate new drag from Lupkes(2015) equations 860 CALL Cdn10_Lupkes2015( ptsui, Cd_ice, Ch_ice )860 CALL Cdn10_Lupkes2015( ptsui, pslp, Cd_ice, Ch_ice ) 861 861 Ce_ice(:,:) = Ch_ice(:,:) ! sensible and latent heat transfer coef. are considered identical 862 862 ENDIF 863 863 864 IF ( iom_use("Cd_ice") ) CALL iom_put("Cd_ice", Cd_ice) ! output value of pure ice-atm. transfer coef.865 IF ( iom_use("Ch_ice") ) CALL iom_put("Ch_ice", Ch_ice) ! output value of pure ice-atm. transfer coef.864 !! IF ( iom_use("Cd_ice") ) CALL iom_put("Cd_ice", Cd_ice) ! output value of pure ice-atm. transfer coef. 865 !! IF ( iom_use("Ch_ice") ) CALL iom_put("Ch_ice", Ch_ice) ! output value of pure ice-atm. transfer coef. 866 866 867 867 ! local scalars ( place there for vector optimisation purposes) 868 ! Computing density of air! Way denser that 1.2 over sea-ice !!! 869 zrhoa (:,:) = rho_air( ptair(:,:), phumi(:,:), pslp(:,:) ) 868 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) 870 869 zcd_dui(:,:) = wndm_ice(:,:) * Cd_ice(:,:) 871 870 … … 877 876 DO jj = 2, jpjm1 878 877 DO ji = fs_2, fs_jpim1 ! vect. opt. 879 putaui(ji,jj) = 0.5_wp * ( zrhoa(ji+1,jj) * zcd_dui(ji+1,jj) &880 & + zrhoa(ji ,jj) * zcd_dui(ji ,jj) ) &878 putaui(ji,jj) = 0.5_wp * ( rhoa(ji+1,jj) * zcd_dui(ji+1,jj) & 879 & + rhoa(ji ,jj) * zcd_dui(ji ,jj) ) & 881 880 & * ( 0.5_wp * ( pwndi(ji+1,jj) + pwndi(ji,jj) ) - rn_vfac * puice(ji,jj) ) 882 pvtaui(ji,jj) = 0.5_wp * ( zrhoa(ji,jj+1) * zcd_dui(ji,jj+1) &883 & + zrhoa(ji,jj ) * zcd_dui(ji,jj ) ) &881 pvtaui(ji,jj) = 0.5_wp * ( rhoa(ji,jj+1) * zcd_dui(ji,jj+1) & 882 & + rhoa(ji,jj ) * zcd_dui(ji,jj ) ) & 884 883 & * ( 0.5_wp * ( pwndj(ji,jj+1) + pwndj(ji,jj) ) - rn_vfac * pvice(ji,jj) ) 885 884 END DO … … 898 897 pevpi (ji,jj) = wndm_ice(ji,jj) * Ce_ice(ji,jj) 899 898 zootm_su = zztmp2 / ptsui(ji,jj) ! ptsui is in K (it can't be zero ??) 900 pssqi (ji,jj) = zztmp1 * EXP( zootm_su ) / zrhoa(ji,jj)899 pssqi (ji,jj) = zztmp1 * EXP( zootm_su ) / rhoa(ji,jj) 901 900 END DO 902 901 END DO … … 934 933 REAL(wp) :: zst3 ! local variable 935 934 REAL(wp) :: zcoef_dqlw, zcoef_dqla ! - - 936 REAL(wp) :: zztmp, zztmp2, z1_rLsub 935 REAL(wp) :: zztmp, zztmp2, z1_rLsub ! - - 937 936 REAL(wp) :: zfr1, zfr2 ! local variables 938 937 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_st ! inverse of surface temperature … … 942 941 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z_dqsb ! sensible heat sensitivity over ice 943 942 REAL(wp), DIMENSION(jpi,jpj) :: zevap, zsnw ! evaporation and snw distribution after wind blowing (SI3) 944 REAL(wp), DIMENSION(jpi,jpj) :: zrhoa 945 !!--------------------------------------------------------------------- 946 ! 947 zcoef_dqlw = 4.0 * 0.95 * Stef ! local scalars 948 zcoef_dqla = -Ls * 11637800. * (-5897.8) 949 ! 950 zrhoa(:,:) = rho_air( ptair(:,:), phumi(:,:), pslp(:,:) ) 943 REAL(wp), DIMENSION(jpi,jpj) :: zqair ! specific humidity of air at z=rn_zqt [kg/kg] !LB 944 !!--------------------------------------------------------------------- 945 ! 946 zcoef_dqlw = 4._wp * 0.95_wp * stefan ! local scalars 947 zcoef_dqla = -rLsub * 11637800._wp * (-5897.8_wp) !LB: BAD! 948 ! 949 SELECT CASE( nhumi ) 950 CASE( np_humi_sph ) 951 zqair(:,:) = phumi(:,:) ! what we read in file is already a spec. humidity! 952 CASE( np_humi_dpt ) 953 zqair(:,:) = q_sat( phumi(:,:), pslp ) 954 CASE( np_humi_rlh ) 955 zqair(:,:) = q_air_rh( 0.01_wp*phumi(:,:), ptair(:,:), pslp(:,:) ) !LB: 0.01 => RH is % percent in file 956 END SELECT 951 957 ! 952 958 zztmp = 1. / ( 1. - albo ) 953 WHERE( ptsu(:,:,:) /= 0._wp ) ; z1_st(:,:,:) = 1._wp / ptsu(:,:,:) 954 ELSEWHERE ; z1_st(:,:,:) = 0._wp 959 WHERE( ptsu(:,:,:) /= 0._wp ) 960 z1_st(:,:,:) = 1._wp / ptsu(:,:,:) 961 ELSEWHERE 962 z1_st(:,:,:) = 0._wp 955 963 END WHERE 956 964 ! ! ========================== ! … … 966 974 qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj) 967 975 ! Long Wave (lw) 968 z_qlw(ji,jj,jl) = 0.95 * ( pqlw(ji,jj) - Stef* ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)976 z_qlw(ji,jj,jl) = 0.95 * ( pqlw(ji,jj) - stefan * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1) 969 977 ! lw sensitivity 970 978 z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3 … … 976 984 ! ... turbulent heat fluxes with Ch_ice recalculated in blk_ice_1 977 985 ! Sensible Heat 978 z_qsb(ji,jj,jl) = zrhoa(ji,jj) * cpa* Ch_ice(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - ptair(ji,jj))986 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)) 979 987 ! Latent Heat 980 988 zztmp2 = EXP( -5897.8 * z1_st(ji,jj,jl) ) 981 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls* Ce_ice(ji,jj) * wndm_ice(ji,jj) * &982 & ( 11637800. * zztmp2 / zrhoa(ji,jj) - phumi(ji,jj) ) )989 qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa(ji,jj) * rLsub * Ce_ice(ji,jj) * wndm_ice(ji,jj) * & 990 & ( 11637800. * zztmp2 / rhoa(ji,jj) - zqair(ji,jj) ) ) 983 991 ! Latent heat sensitivity for ice (Dqla/Dt) 984 992 IF( qla_ice(ji,jj,jl) > 0._wp ) THEN … … 990 998 991 999 ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice) 992 z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa* Ch_ice(ji,jj) * wndm_ice(ji,jj)993 1000 z_dqsb(ji,jj,jl) = rhoa(ji,jj) * rCp_air * Ch_ice(ji,jj) * wndm_ice(ji,jj) 1001 994 1002 ! ----------------------------! 995 1003 ! III Total FLUXES ! … … 1042 1050 DO jl = 1, jpl 1043 1051 qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * rcpi * tmask(:,:,1) ) 1044 ! ! But we do not have Tice => consider it at 0degC => evap=0 1052 ! ! But we do not have Tice => consider it at 0degC => evap=0 1045 1053 END DO 1046 1054 … … 1049 1057 zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) ! zfr2 such that zfr1 + zfr2 to equal 1 1050 1058 ! 1051 WHERE ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) < 0.1_wp ) ! linear decrease from hi=0 to 10cm 1059 WHERE ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) < 0.1_wp ) ! linear decrease from hi=0 to 10cm 1052 1060 qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) ) 1053 1061 ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp ) ! constant (zfr1) when hi>10cm 1054 1062 qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1 1055 1063 ELSEWHERE ! zero when hs>0 1056 qtr_ice_top(:,:,:) = 0._wp 1064 qtr_ice_top(:,:,:) = 0._wp 1057 1065 END WHERE 1058 1066 ! … … 1064 1072 CALL prt_ctl(tab3d_1=ptsu , clinfo1=' blk_ice: ptsu : ', tab3d_2=qns_ice , clinfo2=' qns_ice : ', kdim=jpl) 1065 1073 CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice: tprecip : ', tab2d_2=sprecip , clinfo2=' sprecip : ') 1066 END 1074 ENDIF 1067 1075 ! 1068 1076 END SUBROUTINE blk_ice_2 1069 1077 1070 1078 1071 1079 SUBROUTINE blk_ice_qcn( ld_virtual_itd, ptsu, ptb, phs, phi ) … … 1076 1084 !! to force sea ice / snow thermodynamics 1077 1085 !! in the case conduction flux is emulated 1078 !! 1086 !! 1079 1087 !! ** Method : compute surface energy balance assuming neglecting heat storage 1080 1088 !! following the 0-layer Semtner (1976) approach … … 1101 1109 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zgfac ! enhanced conduction factor 1102 1110 !!--------------------------------------------------------------------- 1103 1111 1104 1112 ! -------------------------------------! 1105 1113 ! I Enhanced conduction factor ! … … 1109 1117 ! 1110 1118 zgfac(:,:,:) = 1._wp 1111 1119 1112 1120 IF( ld_virtual_itd ) THEN 1113 1121 ! … … 1115 1123 zfac2 = EXP(1._wp) * 0.5_wp * zepsilon 1116 1124 zfac3 = 2._wp / zepsilon 1117 ! 1118 DO jl = 1, jpl 1125 ! 1126 DO jl = 1, jpl 1119 1127 DO jj = 1 , jpj 1120 1128 DO ji = 1, jpi … … 1124 1132 END DO 1125 1133 END DO 1126 ! 1127 ENDIF 1128 1134 ! 1135 ENDIF 1136 1129 1137 ! -------------------------------------------------------------! 1130 1138 ! II Surface temperature and conduction flux ! … … 1136 1144 DO jj = 1 , jpj 1137 1145 DO ji = 1, jpi 1138 ! 1146 ! 1139 1147 zkeff_h = zfac * zgfac(ji,jj,jl) / & ! Effective conductivity of the snow-ice system divided by thickness 1140 1148 & ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) ) … … 1153 1161 qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) 1154 1162 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) ) & 1155 1163 & * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) ) 1156 1164 1157 1165 ! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- ! 1158 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) 1166 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) 1159 1167 1160 1168 END DO 1161 1169 END DO 1162 1170 ! 1163 END DO 1164 ! 1171 END DO 1172 ! 1165 1173 END SUBROUTINE blk_ice_qcn 1166 1174 1167 1175 1168 1176 SUBROUTINE Cdn10_Lupkes2012( pcd ) … … 1170 1178 !! *** ROUTINE Cdn10_Lupkes2012 *** 1171 1179 !! 1172 !! ** Purpose : Recompute the neutral air-ice drag referenced at 10m 1180 !! ** Purpose : Recompute the neutral air-ice drag referenced at 10m 1173 1181 !! to make it dependent on edges at leads, melt ponds and flows. 1174 1182 !! After some approximations, this can be resumed to a dependency 1175 1183 !! on ice concentration. 1176 !! 1184 !! 1177 1185 !! ** Method : The parameterization is taken from Lupkes et al. (2012) eq.(50) 1178 1186 !! with the highest level of approximation: level4, eq.(59) … … 1186 1194 !! 1187 1195 !! This new drag has a parabolic shape (as a function of A) starting at 1188 !! Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5 1196 !! Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5 1189 1197 !! and going down to Cdi(say 1.4e-3) for A=1 1190 1198 !! … … 1211 1219 1212 1220 ! ice-atm drag 1213 pcd(:,:) = r cdice + & ! pure ice drag1221 pcd(:,:) = rCd_ice + & ! pure ice drag 1214 1222 & zCe * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp) ! change due to sea-ice morphology 1215 1223 1216 1224 END SUBROUTINE Cdn10_Lupkes2012 1217 1225 1218 1226 1219 SUBROUTINE Cdn10_Lupkes2015( ptm_su, p cd, pch )1227 SUBROUTINE Cdn10_Lupkes2015( ptm_su, pslp, pcd, pch ) 1220 1228 !!---------------------------------------------------------------------- 1221 1229 !! *** ROUTINE Cdn10_Lupkes2015 *** 1222 1230 !! 1223 1231 !! ** pUrpose : Alternative turbulent transfert coefficients formulation 1224 !! between sea-ice and atmosphere with distinct momentum 1225 !! and heat coefficients depending on sea-ice concentration 1232 !! between sea-ice and atmosphere with distinct momentum 1233 !! and heat coefficients depending on sea-ice concentration 1226 1234 !! and atmospheric stability (no meltponds effect for now). 1227 !! 1235 !! 1228 1236 !! ** Method : The parameterization is adapted from Lupkes et al. (2015) 1229 1237 !! and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme, 1230 1238 !! it considers specific skin and form drags (Andreas et al. 2010) 1231 !! to compute neutral transfert coefficients for both heat and 1239 !! to compute neutral transfert coefficients for both heat and 1232 1240 !! momemtum fluxes. Atmospheric stability effect on transfert 1233 1241 !! coefficient is also taken into account following Louis (1979). … … 1239 1247 ! 1240 1248 REAL(wp), DIMENSION(:,:), INTENT(in ) :: ptm_su ! sea-ice surface temperature [K] 1249 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pslp ! sea-level pressure [Pa] 1241 1250 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pcd ! momentum transfert coefficient 1242 1251 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pch ! heat transfert coefficient … … 1267 1276 REAL(wp) :: z0w, z0i, zfmi, zfmw, zfhi, zfhw 1268 1277 REAL(wp) :: zCdn_form_tmp 1269 !!---------------------------------------------------------------------- -----1270 ! 1278 !!---------------------------------------------------------------------- 1279 1271 1280 ! Momentum Neutral Transfert Coefficients (should be a constant) 1272 1281 zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2 ! Eq. 40 … … 1277 1286 ! Heat Neutral Transfert Coefficients 1278 1287 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 1279 1288 1280 1289 ! Atmospheric and Surface Variables 1281 1290 zst(:,:) = sst_m(:,:) + rt0 ! convert SST from Celcius to Kelvin 1282 zqo_sat(:,:) = 0.98_wp * q_sat( zst(:,:) , sf(jp_slp)%fnow(:,:,1) )! saturation humidity over ocean [kg/kg]1283 zqi_sat(:,:) = 0.98_wp * q_sat( ptm_su(:,:), sf(jp_slp)%fnow(:,:,1) )! saturation humidity over ice [kg/kg]1284 ! 1285 DO jj = 2, jpjm1 ! reduced loop is necessary for reproduc tibility1291 zqo_sat(:,:) = rdct_qsat_salt * q_sat( zst(:,:) , pslp(:,:) ) ! saturation humidity over ocean [kg/kg] 1292 zqi_sat(:,:) = q_sat( ptm_su(:,:), pslp(:,:) ) ! saturation humidity over ice [kg/kg] 1293 ! 1294 DO jj = 2, jpjm1 ! reduced loop is necessary for reproducibility 1286 1295 DO ji = fs_2, fs_jpim1 1287 1296 ! Virtual potential temperature [K] … … 1289 1298 zthetav_is = ptm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) ) ! ocean ice 1290 1299 zthetav_zu = t_zu (ji,jj) * ( 1._wp + rctv0 * q_zu(ji,jj) ) ! at zu 1291 1300 1292 1301 ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead) 1293 1302 zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj) )**2 ! over ocean 1294 1303 zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2 ! over ice 1295 1304 1296 1305 ! Momentum and Heat Neutral Transfert Coefficients 1297 1306 zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta ! Eq. 40 1298 zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) ) ! Eq. 53 1299 1300 ! Momentum and Heat Stability functions ( possibility to use psi_m_ecmwf instead)1307 zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) ) ! Eq. 53 1308 1309 ! Momentum and Heat Stability functions (!!GS: possibility to use psi_m_ecmwf instead ?) 1301 1310 z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water 1302 1311 z0i = z0_skin_ice ! over ice … … 1309 1318 zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) ) ! Eq. 28 1310 1319 ENDIF 1311 1320 1312 1321 IF( zrib_i <= 0._wp ) THEN 1313 1322 zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp))) ! Eq. 9 … … 1317 1326 zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) ) ! Eq. 27 1318 1327 ENDIF 1319 1328 1320 1329 ! Momentum Transfert Coefficients (Eq. 38) 1321 1330 pcd(ji,jj) = zCdn_skin_ice * zfmi + & 1322 1331 & 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) ) 1323 1332 1324 1333 ! Heat Transfert Coefficients (Eq. 49) 1325 1334 pch(ji,jj) = zChn_skin_ice * zfhi + &
Note: See TracChangeset
for help on using the changeset viewer.
