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