- Timestamp:
- 2010-04-30T17:49:04+02:00 (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/DEV_r1837_mass_heat_salt_fluxes/NEMO/LIM_SRC_2/limthd_zdf_2.F90
r1756 r1855 4 4 !! thermodynamic growth and decay of the ice 5 5 !!====================================================================== 6 !! History : 1.0 ! 01-04 (LIM) Original code7 !! 2.0 ! 02-08 (C. Ethe, G. Madec) F906 !! History : 1.0 ! 2001-04 (LIM) Original code 7 !! 2.0 ! 2002-08 (C. Ethe, G. Madec) F90 8 8 !!---------------------------------------------------------------------- 9 9 #if defined key_lim2 … … 11 11 !! 'key_lim2' LIM 2.0 sea-ice model 12 12 !!---------------------------------------------------------------------- 13 !!----------------------------------------------------------------------14 13 !! lim_thd_zdf_2 : vertical accr./abl. and lateral ablation of sea ice 15 14 !!---------------------------------------------------------------------- 16 !! * Modules used17 15 USE par_oce ! ocean parameters 18 USE phycst ! ???19 USE thd_ice_220 USE ice_221 USE limistate_222 USE in_out_manager16 USE phycst ! physical constants 17 USE ice_2 ! LIM 2 sea-ice variables 18 USE limistate_2 ! LIM 2 initial state 19 USE in_out_manager ! I/O manager 20 !! 23 21 USE cpl_oasis3, ONLY : lk_cpl 24 22 … … 26 24 PRIVATE 27 25 28 PUBLIC lim_thd_zdf_2 ! called by lim_thd_2 29 30 REAL(wp) :: epsi20 = 1.e-20 , & ! constant values 31 & epsi13 = 1.e-13 , & 32 & zzero = 0.e0 , & 33 & zone = 1.e0 26 PUBLIC lim_thd_zdf_2 ! called by lim_thd_2 27 28 REAL(wp) :: epsi20 = 1.e-20 ! constant values 29 REAL(wp) :: epsi13 = 1.e-13 ! 30 REAL(wp) :: rzero = 0.e0 ! 31 REAL(wp) :: rone = 1.e0 ! 32 34 33 !!---------------------------------------------------------------------- 35 !! LIM 2.0, UCL-LOCEAN-IPSL (2005)34 !! NEMO/LIM 3.3, UCL-LOCEAN-IPSL (2010) 36 35 !! $Id$ 37 36 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 73 72 !! 74 73 INTEGER :: ji ! dummy loop indices 75 REAL(wp), DIMENSION(jpij,2) :: zqcmlt ! energy due to surface( /1 ) and bottom melting( /2 ) 76 REAL(wp), DIMENSION(jpij) :: & 77 ztsmlt & ! snow/ice surface melting temperature 78 ,ztbif & ! int. temp. at the mid-point of the 1st layer of the snow/ice sys. 79 ,zksn & ! effective conductivity of snow 80 ,zkic & ! effective conductivity of ice 81 ,zksndh & ! thermal cond. at the mid-point of the 1st layer of the snow/ice sys. 82 , zfcsu & ! conductive heat flux at the surface of the snow/ice system 83 , zfcsudt & ! = zfcsu * dt 84 , zi0 & ! frac. of the net SW rad. which is not absorbed at the surface 85 , z1mi0 & ! fraction of the net SW radiation absorbed at the surface 86 , zqmax & ! maximum energy stored in brine pockets 87 , zrcpdt & ! h_su*rho_su*cp_su/dt(h_su being the thick. of surf. layer) 88 , zts_old & ! previous surface temperature 89 , zidsn , z1midsn , zidsnic ! tempory variables 90 REAL(wp), DIMENSION(jpij) :: & 91 zfnet & ! net heat flux at the top surface( incl. conductive heat flux) 92 , zsprecip & ! snow accumulation 93 , zhsnw_old & ! previous snow thickness 94 , zdhictop & ! change in ice thickness due to top surf ablation/accretion 95 , zdhicbot & ! change in ice thickness due to bottom surf abl/acc 96 , zqsup & ! energy transmitted to ocean (coming from top surf abl/acc) 97 , zqocea & ! energy transmitted to ocean (coming from bottom sur abl/acc) 98 , zfrl_old & ! previous sea/ice fraction 99 , zfrld_1d & ! new sea/ice fraction 100 , zep ! internal temperature of the 2nd layer of the snow/ice system 101 REAL(wp), DIMENSION(3) :: & 102 zplediag & ! principle diagonal, subdiag. and supdiag. of the 103 , zsubdiag & ! tri-diagonal matrix coming from the computation 104 , zsupdiag & ! of the temperatures inside the snow-ice system 105 , zsmbr ! second member 106 REAL(wp) :: & 107 zhsu & ! thickness of surface layer 108 , zhe & ! effective thickness for compu. of equ. thermal conductivity 109 , zheshth & ! = zhe / thth 110 , zghe & ! correction factor of the thermal conductivity 111 , zumsb & ! parameter for numerical method to solve heat-diffusion eq. 112 , zkhsn & ! conductivity at the snow layer 113 , zkhic & ! conductivity at the ice layers 114 , zkint & ! equivalent conductivity at the snow-ice interface 115 , zkhsnint & ! = zkint*dt / (hsn*rhosn*cpsn) 116 , zkhicint & ! = 2*zkint*dt / (hic*rhoic*cpic) 117 , zpiv1 , zpiv2 & ! tempory scalars used to solve the tri-diagonal system 118 , zb2 , zd2 , zb3 , zd3 & 119 , ztint ! equivalent temperature at the snow-ice interface 120 REAL(wp) :: & 121 zexp & ! exponential function of the ice thickness 122 , zfsab & ! part of solar radiation stored in brine pockets 123 , zfts & ! value of energy balance function when the temp. equal surf. temp. 124 , zdfts & ! value of derivative of ztfs when the temp. equal surf. temp. 125 , zdts & ! surface temperature increment 126 , zqsnw_mlt & ! energy needed to melt snow 127 , zdhsmlt & ! change in snow thickness due to melt 128 , zhsn & ! snow thickness (previous+accumulation-melt) 129 , zqsn_mlt_rem & ! remaining heat coming from snow melting 130 , zqice_top_mlt & ! energy used to melt ice at top surface 131 , zdhssub & ! change in snow thick. due to sublimation or evaporation 132 , zdhisub & ! change in ice thick. due to sublimation or evaporation 133 , zdhsn & ! snow ice thickness increment 134 , zdtsn & ! snow internal temp. increment 135 , zdtic & ! ice internal temp. increment 136 , zqnes ! conductive energy due to ice melting in the first ice layer 137 REAL(wp) :: & 138 ztbot & ! temperature at the bottom surface 139 , zfcbot & ! conductive heat flux at bottom surface 140 , zqice_bot & ! energy used for bottom melting/growing 141 , zqice_bot_mlt & ! energy used for bottom melting 142 , zqstbif_bot & ! part of energy stored in brine pockets used for bottom melting 143 , zqstbif_old & ! tempory var. for zqstbif_bot 144 , zdhicmlt & ! change in ice thickness due to bottom melting 145 , zdhicm & ! change in ice thickness var. 146 , zdhsnm & ! change in snow thickness var. 147 , zhsnfi & ! snow thickness var. 148 , zc1, zpc1, zc2, zpc2, zp1, zp2 & ! tempory variables 149 , ztb2, ztb3 150 REAL(wp) :: & 151 zdrmh & ! change in snow/ice thick. after snow-ice formation 152 , zhicnew & ! new ice thickness 153 , zhsnnew & ! new snow thickness 154 , zquot , ztneq & ! tempory temp. variables 155 , zqice, zqicetot & ! total heat inside the snow/ice system 156 , zdfrl & ! change in ice concentration 157 , zdvsnvol & ! change in snow volume 158 , zdrfrl1, zdrfrl2 & ! tempory scalars 159 , zihsn, zidhb, zihic, zihe, zihq, ziexp, ziqf, zihnf, zibmlt, ziqr, zihgnew, zind 160 !!---------------------------------------------------------------------- 161 162 !----------------------------------------------------------------------- 163 ! 1. Boundaries conditions for snow/ice system internal temperature 164 ! - If tbif_1d(ji,1) > rt0_snow, tbif_1d(ji,1) = rt0_snow 165 ! - If tbif_1d(ji,2/3) > rt0_ice, tbif_1d(ji,2/3) = rt0_ice 166 ! Computation of energies due to surface and bottom melting 167 !----------------------------------------------------------------------- 168 169 DO ji = kideb , kiut 170 zihsn = MAX( zzero , SIGN( zone , hsndif - h_snow_1d(ji) ) ) 171 zihic = MAX( zzero , SIGN( zone , hicdif - h_ice_1d(ji) ) ) 172 !--computation of energy due to surface melting 173 zqcmlt(ji,1) = ( MAX ( zzero , & 174 & rcpsn * h_snow_1d(ji) * ( tbif_1d(ji,1) - rt0_snow ) ) ) * ( 1.0 - zihsn ) 175 !--computation of energy due to bottom melting 176 zqcmlt(ji,2) = ( MAX( zzero , & 177 & rcpic * ( tbif_1d(ji,2) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) ) & 178 & + MAX( zzero , & 179 & rcpic * ( tbif_1d(ji,3) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) ) & 180 & ) * ( 1.0 - zihic ) 181 !--limitation of snow/ice system internal temperature 182 tbif_1d(ji,1) = MIN( rt0_snow, tbif_1d(ji,1) ) 183 tbif_1d(ji,2) = MIN( rt0_ice , tbif_1d(ji,2) ) 184 tbif_1d(ji,3) = MIN( rt0_ice , tbif_1d(ji,3) ) 185 END DO 186 187 !------------------------------------------- 188 ! 2. Calculate some intermediate variables. 189 !------------------------------------------- 190 191 ! initialisation of the thickness of surface layer 192 zhsu = hnzst 193 194 DO ji = kideb , kiut 195 zind = MAX( zzero , SIGN( zone , zhsu - h_snow_1d(ji) ) ) 196 zihsn = MAX( zzero , SIGN( zone , hsndif - h_snow_1d(ji) ) ) 197 zihsn = MAX( zihsn , zind ) 198 zihic = MAX( zzero , sign( zone , hicdif - h_ice_1d(ji) ) ) 199 ! 2.1. Computation of surface melting temperature 200 !---------------------------------------------------- 201 zind = MAX( zzero , SIGN( zone , -h_snow_1d(ji) ) ) 202 ztsmlt(ji) = ( 1.0 - zind ) * rt0_snow + zind * rt0_ice 203 ! 204 ! 2.2. Effective conductivity of snow and ice 205 !----------------------------------------------- 206 207 !---computation of the correction factor on the thermal conductivity 208 !-- (Morales Maqueda, 1995 ; Fichefet and Morales Maqueda, 1997) 209 zhe = ( rcdsn / ( rcdsn + rcdic ) ) * h_ice_1d(ji) & 210 & + ( rcdic / ( rcdsn + rcdic ) ) * h_snow_1d(ji) 211 zihe = MAX( zzero , SIGN( zone , 2.0 * zhe - thth ) ) 212 zheshth = zhe / thth 213 zghe = ( 1.0 - zihe ) * zheshth * ( 2.0 - zheshth ) & 214 & + zihe * 0.5 * ( 1.5 + LOG( 2.0 * zheshth ) ) 215 216 !---effective conductivities 217 zksn(ji) = zghe * rcdsn 218 zkic(ji) = zghe * rcdic 219 220 ! 221 ! 2.3. Computation of the conductive heat flux from the snow/ice 222 ! system interior toward the top surface 223 !------------------------------------------------------------------ 224 225 !---Thermal conductivity at the mid-point of the first snow/ice system layer 226 zksndh(ji) = ( ( 1.0 - zihsn ) * 2.0 * zksn(ji) + zihsn * 4.0 * zkic(ji) ) & 227 & / ( ( 1.0 - zihsn ) * h_snow_1d(ji) & 228 & + zihsn * ( ( 1.0 + 3.0 * zihic ) * h_ice_1d(ji) & 229 & + 4.0 * zkic(ji)/zksn(ji) * h_snow_1d(ji) ) ) 230 231 !---internal temperature at the mid-point of the first snow/ice system layer 232 ztbif(ji) = ( 1.0 - zihsn ) * tbif_1d(ji,1) & 233 & + zihsn * ( ( 1.0 - zihic ) * tbif_1d(ji,2) & 234 & + zihic * tfu_1d(ji) ) 235 !---conductive heat flux 236 zfcsu(ji) = zksndh(ji) * ( ztbif(ji) - sist_1d(ji) ) 237 238 END DO 239 240 !-------------------------------------------------------------------- 241 ! 3. Calculate : 242 ! - fstbif_1d, part of solar radiation absorbing inside the ice 243 ! assuming an exponential absorption (Grenfell and Maykut, 1977) 244 ! - zqmax, maximum energy stored in brine pockets 245 ! - qstbif_1d, total energy stored in brine pockets (updating) 246 !------------------------------------------------------------------- 247 248 DO ji = kideb , kiut 249 zihsn = MAX( zzero , SIGN (zone , -h_snow_1d(ji) ) ) 250 zihic = MAX( zzero , 1.0 - ( h_ice_1d(ji) / zhsu ) ) 251 zind = MAX( zzero , SIGN (zone , hicdif - h_ice_1d(ji) ) ) 252 !--Computation of the fraction of the net shortwave radiation which 253 !--penetrates inside the ice cover ( See Forcat) 254 zi0(ji) = zihsn * ( fr1_i0_1d(ji) + zihic * fr2_i0_1d(ji) ) 255 zexp = MIN( zone , EXP( -1.5 * ( h_ice_1d(ji) - zhsu ) ) ) 256 fstbif_1d(ji) = zi0(ji) * qsr_ice_1d(ji) * zexp 257 !--Computation of maximum energy stored in brine pockets zqmax and update 258 !--the total energy stored in brine pockets, if less than zqmax 259 zqmax(ji) = MAX( zzero , 0.5 * xlic * ( h_ice_1d(ji) - hicmin ) ) 260 zfsab = zi0(ji) * qsr_ice_1d(ji) * ( 1.0 - zexp ) 261 zihq = ( 1.0 - zind ) * MAX(zzero, SIGN( zone , qstbif_1d(ji) - zqmax(ji) ) ) & 262 & + zind * zone 263 qstbif_1d(ji) = ( qstbif_1d(ji) + ( 1.0 - zihq ) * zfsab * rdt_ice ) * swiqst 264 !--fraction of shortwave radiation absorbed at surface 265 ziexp = zihq * zexp + ( 1.0 - zihq ) * ( swiqst + ( 1.0 - swiqst ) * zexp ) 266 z1mi0(ji) = 1.0 - zi0(ji) * ziexp 267 END DO 268 269 !-------------------------------------------------------------------------------- 270 ! 4. Computation of the surface temperature : determined by considering the 271 ! budget of a thin layer of thick. zhsu at the top surface (H. Grenier, 1995) 272 ! and based on a surface energy balance : 273 ! hsu * rcp * dT/dt = Fsr + Fnsr(T) + Fcs(T), 274 ! where - Fsr is the net absorbed solar radiation, 275 ! - Fnsr is the total non solar radiation (incoming and outgoing long-wave, 276 ! sensible and latent heat fluxes) 277 ! - Fcs the conductive heat flux at the top of surface 278 !------------------------------------------------------------------------------ 279 280 ! 4.1. Computation of intermediate values 281 !--------------------------------------------- 282 DO ji = kideb, kiut 283 zrcpdt(ji) = ( rcpsn * MIN( h_snow_1d(ji) , zhsu ) & 284 & + rcpic * MAX( zhsu - h_snow_1d(ji) , zzero ) ) / rdt_ice 285 zts_old(ji) = sist_1d(ji) 286 END DO 287 288 ! 4.2. Computation of surface temperature by expanding the eq. of energy balance 289 ! with Ts = Tp + DT. One obtain , F(Tp) + DT * DF(Tp) = 0 290 ! where - F(Tp) = Fsr + Fnsr(Tp) + Fcs(Tp) 291 ! - DF(Tp)= (dFnsr(Tp)/dT) + (dFcs(Tp)/dT) - hsu*rcp/dt 292 !--------------------------------------------------------------------------------- 293 294 DO ji = kideb, kiut 295 !---computation of the derivative of energy balance function 296 zdfts = zksndh(ji) & ! contribution of the conductive heat flux 297 & + zrcpdt(ji) & ! contribution of hsu * rcp / dt 298 & - dqns_ice_1d (ji) ! contribution of the total non solar radiation 299 !---computation of the energy balance function 300 zfts = - z1mi0 (ji) * qsr_ice_1d(ji) & ! net absorbed solar radiation 301 & - qns_ice_1d(ji) & ! total non solar radiation 302 & - zfcsu (ji) ! conductive heat flux from the surface 303 !---computation of surface temperature increment 304 zdts = -zfts / zdfts 305 !---computation of the new surface temperature 306 sist_1d(ji) = sist_1d(ji) + zdts 307 END DO 308 309 !---------------------------------------------------------------------------- 310 ! 5. Boundary condition at the top surface 311 !-- IF Tsb < Tmelt, Fnet = Fcs (the net heat flux equal the conductive heat flux) 312 ! Otherwise Tsb = Tmelt and Qnet(Tmelt) > 0 313 ! Fnet(Tmelt) is therefore the net surface flux needed for melting 314 !---------------------------------------------------------------------------- 74 REAL(wp), DIMENSION(jpij,2) :: zqcmlt ! energy due to surface( /1 ) and bottom melting( /2 ) 75 REAL(wp), DIMENSION(jpij) :: ztsmlt ! snow/ice surface melting temperature 76 REAL(wp), DIMENSION(jpij) :: ztbif ! int. temp. at the mid-point of the 1st layer of the snow/ice sys. 77 REAL(wp), DIMENSION(jpij) :: zksn ! effective conductivity of snow 78 REAL(wp), DIMENSION(jpij) :: zkic ! effective conductivity of ice 79 REAL(wp), DIMENSION(jpij) :: zksndh ! thermal cond. at the mid-point of the 1st layer of the snow/ice sys. 80 REAL(wp), DIMENSION(jpij) :: zfcsu ! conductive heat flux at the surface of the snow/ice system 81 REAL(wp), DIMENSION(jpij) :: zfcsudt ! = zfcsu * dt 82 REAL(wp), DIMENSION(jpij) :: zi0 ! frac. of the net SW rad. which is not absorbed at the surface 83 REAL(wp), DIMENSION(jpij) :: z1mi0 ! fraction of the net SW radiation absorbed at the surface 84 REAL(wp), DIMENSION(jpij) :: zqmax ! maximum energy stored in brine pockets 85 REAL(wp), DIMENSION(jpij) :: zrcpdt ! h_su*rho_su*cp_su/dt(h_su being the thick. of surf. layer) 86 REAL(wp), DIMENSION(jpij) :: zts_old ! previous surface temperature 87 REAL(wp), DIMENSION(jpij) :: zidsn , z1midsn , zidsnic ! tempory variables 88 !! 89 REAL(wp), DIMENSION(jpij) :: zfnet ! net heat flux at the top surface( incl. conductive heat flux) 90 REAL(wp), DIMENSION(jpij) :: zsprecip ! snow accumulation 91 REAL(wp), DIMENSION(jpij) :: zhsnw_old ! previous snow thickness 92 REAL(wp), DIMENSION(jpij) :: zdhictop ! change in ice thickness due to top surf ablation/accretion 93 REAL(wp), DIMENSION(jpij) :: zdhicbot ! change in ice thickness due to bottom surf abl/acc 94 REAL(wp), DIMENSION(jpij) :: zqsup ! energy transmitted to ocean (coming from top surf abl/acc) 95 REAL(wp), DIMENSION(jpij) :: zqocea ! energy transmitted to ocean (coming from bottom sur abl/acc) 96 REAL(wp), DIMENSION(jpij) :: zfrl_old ! previous sea/ice fraction 97 REAL(wp), DIMENSION(jpij) :: zfrld_1d ! new sea/ice fraction 98 REAL(wp), DIMENSION(jpij) :: zep ! internal temperature of the 2nd layer of the snow/ice system 99 !! 100 REAL(wp), DIMENSION(3) :: zplediag ! principle diagonal, subdiag. and supdiag. of the 101 REAL(wp), DIMENSION(3) :: zsubdiag ! tri-diagonal matrix coming from the computation 102 REAL(wp), DIMENSION(3) :: zsupdiag ! of the temperatures inside the snow-ice system 103 REAL(wp), DIMENSION(3) :: zsmbr ! second member 104 REAL(wp) :: zhsu ! thickness of surface layer 105 REAL(wp) :: zhe ! effective thickness for compu. of equ. thermal conductivity 106 REAL(wp) :: zheshth ! = zhe / thth 107 REAL(wp) :: zghe ! correction factor of the thermal conductivity 108 REAL(wp) :: zumsb ! parameter for numerical method to solve heat-diffusion eq. 109 REAL(wp) :: zkhsn ! conductivity at the snow layer 110 REAL(wp) :: zkhic ! conductivity at the ice layers 111 REAL(wp) :: zkint ! equivalent conductivity at the snow-ice interface 112 REAL(wp) :: zkhsnint ! = zkint*dt / (hsn*rhosn*cpsn) 113 REAL(wp) :: zkhicint ! = 2*zkint*dt / (hic*rhoic*cpic) 114 REAL(wp) :: zpiv1 , zpiv2 ! tempory scalars used to solve the tri-diagonal system 115 REAL(wp) :: zb2, zd2, zb3, zd3 116 REAL(wp) :: ztint ! equivalent temperature at the snow-ice interface 117 REAL(wp) :: zexp ! exponential function of the ice thickness 118 REAL(wp) :: zfsab ! part of solar radiation stored in brine pockets 119 REAL(wp) :: zfts ! value of energy balance function when the temp. equal surf. temp. 120 REAL(wp) :: zdfts ! value of derivative of ztfs when the temp. equal surf. temp. 121 REAL(wp) :: zdts ! surface temperature increment 122 REAL(wp) :: zqsnw_mlt ! energy needed to melt snow 123 REAL(wp) :: zdhsmlt ! change in snow thickness due to melt 124 REAL(wp) :: zhsn ! snow thickness (previous+accumulation-melt) 125 REAL(wp) :: zqsn_mlt_rem ! remaining heat coming from snow melting 126 REAL(wp) :: zqice_top_mlt ! energy used to melt ice at top surface 127 REAL(wp) :: zdhssub ! change in snow thick. due to sublimation or evaporation 128 REAL(wp) :: zdhisub ! change in ice thick. due to sublimation or evaporation 129 REAL(wp) :: zdhsn ! snow ice thickness increment 130 REAL(wp) :: zdtsn ! snow internal temp. increment 131 REAL(wp) :: zdtic ! ice internal temp. increment 132 REAL(wp) :: zqnes ! conductive energy due to ice melting in the first ice layer 133 !! 134 REAL(wp) :: ztbot ! temperature at the bottom surface 135 REAL(wp) :: zfcbot ! conductive heat flux at bottom surface 136 REAL(wp) :: zqice_bot ! energy used for bottom melting/growing 137 REAL(wp) :: zqice_bot_mlt ! energy used for bottom melting 138 REAL(wp) :: zqstbif_bot ! part of energy stored in brine pockets used for bottom melting 139 REAL(wp) :: zqstbif_old ! tempory var. for zqstbif_bot 140 REAL(wp) :: zdhicmlt ! change in ice thickness due to bottom melting 141 REAL(wp) :: zdhicm ! change in ice thickness var. 142 REAL(wp) :: zdhsnm ! change in snow thickness var. 143 REAL(wp) :: zhsnfi ! snow thickness var. 144 REAL(wp) :: zc1, zpc1, zc2, zpc2, zp1, zp2 ! temporary scalars 145 REAL(wp) :: ztb2, ztb3 146 !! 147 REAL(wp) :: zdrmh ! change in snow/ice thick. after snow-ice formation 148 REAL(wp) :: zhicnew ! new ice thickness 149 REAL(wp) :: zhsnnew ! new snow thickness 150 REAL(wp) :: zquot , ztneq ! tempory temp. variables 151 REAL(wp) :: zqice, zqicetot & ! total heat inside the snow/ice system 152 REAL(wp) :: zdfrl ! change in ice concentration 153 REAL(wp) :: zdvsnvol ! change in snow volume 154 REAL(wp) :: zdrfrl1, zdrfrl2 ! tempory scalars 155 REAL(wp) :: zihsn, zidhb, zihic, zihe, zihq, ziexp, ziqf, zihnf, zibmlt, ziqr, zihgnew, zind 156 !!---------------------------------------------------------------------- 157 158 !----------------------------------------------------------------------- 159 ! 1. Boundaries conditions for snow/ice system internal temperature 160 ! - If tbif_1d(ji,1) > rt0_snow, tbif_1d(ji,1) = rt0_snow 161 ! - If tbif_1d(ji,2/3) > rt0_ice, tbif_1d(ji,2/3) = rt0_ice 162 ! Computation of energies due to surface and bottom melting 163 !----------------------------------------------------------------------- 164 DO ji = kideb , kiut 165 zihsn = MAX( zzero , SIGN( zone , hsndif - h_snow_1d(ji) ) ) 166 zihic = MAX( zzero , SIGN( zone , hicdif - h_ice_1d(ji) ) ) 167 !--computation of energy due to surface melting 168 zqcmlt(ji,1) = ( MAX ( zzero , & 169 & rcpsn * h_snow_1d(ji) * ( tbif_1d(ji,1) - rt0_snow ) ) ) * ( 1.0 - zihsn ) 170 !--computation of energy due to bottom melting 171 zqcmlt(ji,2) = ( MAX( zzero , & 172 & rcpic * ( tbif_1d(ji,2) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) ) & 173 & + MAX( zzero , & 174 & rcpic * ( tbif_1d(ji,3) - rt0_ice ) * ( h_ice_1d(ji) / 2. ) ) & 175 & ) * ( 1.0 - zihic ) 176 !--limitation of snow/ice system internal temperature 177 tbif_1d(ji,1) = MIN( rt0_snow, tbif_1d(ji,1) ) 178 tbif_1d(ji,2) = MIN( rt0_ice , tbif_1d(ji,2) ) 179 tbif_1d(ji,3) = MIN( rt0_ice , tbif_1d(ji,3) ) 180 END DO 181 182 !------------------------------------------- 183 ! 2. Calculate some intermediate variables. 184 !------------------------------------------- 185 zhsu = hnzst ! initialisation of the thickness of surface layer 186 ! 187 DO ji = kideb , kiut 188 zind = MAX( zzero , SIGN( zone , zhsu - h_snow_1d(ji) ) ) 189 zihsn = MAX( zzero , SIGN( zone , hsndif - h_snow_1d(ji) ) ) 190 zihsn = MAX( zihsn , zind ) 191 zihic = MAX( zzero , sign( zone , hicdif - h_ice_1d (ji) ) ) 192 ! 193 ! 2.1. Computation of surface melting temperature 194 !---------------------------------------------------- 195 zind = MAX( zzero , SIGN( zone , -h_snow_1d(ji) ) ) 196 ztsmlt(ji) = ( 1.0 - zind ) * rt0_snow + zind * rt0_ice 197 ! 198 ! 2.2. Effective conductivity of snow and ice 199 !----------------------------------------------- 200 201 !---computation of the correction factor on the thermal conductivity 202 !-- (Morales Maqueda, 1995 ; Fichefet and Morales Maqueda, 1997) 203 zhe = ( rcdsn / ( rcdsn + rcdic ) ) * h_ice_1d(ji) & 204 & + ( rcdic / ( rcdsn + rcdic ) ) * h_snow_1d(ji) 205 zihe = MAX( zzero , SIGN( zone , 2.0 * zhe - thth ) ) 206 zheshth = zhe / thth 207 zghe = ( 1.0 - zihe ) * zheshth * ( 2.0 - zheshth ) & 208 & + zihe * 0.5 * ( 1.5 + LOG( 2.0 * zheshth ) ) 209 210 !---effective conductivities 211 zksn(ji) = zghe * rcdsn 212 zkic(ji) = zghe * rcdic 213 ! 214 ! 2.3. Computation of the conductive heat flux from the snow/ice 215 ! system interior toward the top surface 216 !------------------------------------------------------------------ 217 218 !---Thermal conductivity at the mid-point of the first snow/ice system layer 219 zksndh(ji) = ( ( 1.0 - zihsn ) * 2.0 * zksn(ji) + zihsn * 4.0 * zkic(ji) ) & 220 & / ( ( 1.0 - zihsn ) * h_snow_1d(ji) & 221 & + zihsn * ( ( 1.0 + 3.0 * zihic ) * h_ice_1d(ji) & 222 & + 4.0 * zkic(ji)/zksn(ji) * h_snow_1d(ji) ) ) 223 224 !---internal temperature at the mid-point of the first snow/ice system layer 225 ztbif(ji) = ( 1.0 - zihsn ) * tbif_1d(ji,1) & 226 & + zihsn * ( ( 1.0 - zihic ) * tbif_1d(ji,2) & 227 & + zihic * tfu_1d(ji) ) 228 !---conductive heat flux 229 zfcsu(ji) = zksndh(ji) * ( ztbif(ji) - sist_1d(ji) ) 230 231 END DO 232 233 !-------------------------------------------------------------------- 234 ! 3. Calculate : 235 ! - fstbif_1d, part of solar radiation absorbing inside the ice 236 ! assuming an exponential absorption (Grenfell and Maykut, 1977) 237 ! - zqmax, maximum energy stored in brine pockets 238 ! - qstbif_1d, total energy stored in brine pockets (updating) 239 !------------------------------------------------------------------- 240 241 DO ji = kideb , kiut 242 zihsn = MAX( zzero , SIGN (zone , -h_snow_1d(ji) ) ) 243 zihic = MAX( zzero , 1.0 - ( h_ice_1d(ji) / zhsu ) ) 244 zind = MAX( zzero , SIGN (zone , hicdif - h_ice_1d(ji) ) ) 245 !--Computation of the fraction of the net shortwave radiation which 246 !--penetrates inside the ice cover ( See Forcat) 247 zi0(ji) = zihsn * ( fr1_i0_1d(ji) + zihic * fr2_i0_1d(ji) ) 248 zexp = MIN( zone , EXP( -1.5 * ( h_ice_1d(ji) - zhsu ) ) ) 249 fstbif_1d(ji) = zi0(ji) * qsr_ice_1d(ji) * zexp 250 !--Computation of maximum energy stored in brine pockets zqmax and update 251 !--the total energy stored in brine pockets, if less than zqmax 252 zqmax(ji) = MAX( zzero , 0.5 * xlic * ( h_ice_1d(ji) - hicmin ) ) 253 zfsab = zi0(ji) * qsr_ice_1d(ji) * ( 1.0 - zexp ) 254 zihq = ( 1.0 - zind ) * MAX(zzero, SIGN( zone , qstbif_1d(ji) - zqmax(ji) ) ) & 255 & + zind * zone 256 qstbif_1d(ji) = ( qstbif_1d(ji) + ( 1.0 - zihq ) * zfsab * rdt_ice ) * swiqst 257 !--fraction of shortwave radiation absorbed at surface 258 ziexp = zihq * zexp + ( 1.0 - zihq ) * ( swiqst + ( 1.0 - swiqst ) * zexp ) 259 z1mi0(ji) = 1.0 - zi0(ji) * ziexp 260 END DO 261 262 !-------------------------------------------------------------------------------- 263 ! 4. Computation of the surface temperature : determined by considering the 264 ! budget of a thin layer of thick. zhsu at the top surface (H. Grenier, 1995) 265 ! and based on a surface energy balance : 266 ! hsu * rcp * dT/dt = Fsr + Fnsr(T) + Fcs(T), 267 ! where - Fsr is the net absorbed solar radiation, 268 ! - Fnsr is the total non solar radiation (incoming and outgoing long-wave, 269 ! sensible and latent heat fluxes) 270 ! - Fcs the conductive heat flux at the top of surface 271 !------------------------------------------------------------------------------ 272 273 ! 4.1. Computation of intermediate values 274 !--------------------------------------------- 275 DO ji = kideb, kiut 276 zrcpdt(ji) = ( rcpsn * MIN( h_snow_1d(ji) , zhsu ) & 277 & + rcpic * MAX( zhsu - h_snow_1d(ji) , zzero ) ) / rdt_ice 278 zts_old(ji) = sist_1d(ji) 279 END DO 280 281 ! 4.2. Computation of surface temperature by expanding the eq. of energy balance 282 ! with Ts = Tp + DT. One obtain , F(Tp) + DT * DF(Tp) = 0 283 ! where - F(Tp) = Fsr + Fnsr(Tp) + Fcs(Tp) 284 ! - DF(Tp)= (dFnsr(Tp)/dT) + (dFcs(Tp)/dT) - hsu*rcp/dt 285 !--------------------------------------------------------------------------------- 286 287 DO ji = kideb, kiut 288 !---computation of the derivative of energy balance function 289 zdfts = zksndh(ji) & ! contribution of the conductive heat flux 290 & + zrcpdt(ji) & ! contribution of hsu * rcp / dt 291 & - dqns_ice_1d (ji) ! contribution of the total non solar radiation 292 !---computation of the energy balance function 293 zfts = - z1mi0 (ji) * qsr_ice_1d(ji) & ! net absorbed solar radiation 294 & - qns_ice_1d(ji) & ! total non solar radiation 295 & - zfcsu (ji) ! conductive heat flux from the surface 296 !---computation of surface temperature increment 297 zdts = -zfts / zdfts 298 !---computation of the new surface temperature 299 sist_1d(ji) = sist_1d(ji) + zdts 300 END DO 301 302 !---------------------------------------------------------------------------- 303 ! 5. Boundary condition at the top surface 304 !-- IF Tsb < Tmelt, Fnet = Fcs (the net heat flux equal the conductive heat flux) 305 ! Otherwise Tsb = Tmelt and Qnet(Tmelt) > 0 306 ! Fnet(Tmelt) is therefore the net surface flux needed for melting 307 !---------------------------------------------------------------------------- 315 308 316 309 317 318 319 310 ! 5.1. Limitation of surface temperature and update total non solar fluxes, 311 ! latent heat flux and conductive flux at the top surface 312 !---------------------------------------------------------------------- 320 313 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 314 IF ( .NOT. lk_cpl ) THEN ! duplicate the loop for performances issues 315 DO ji = kideb, kiut 316 sist_1d(ji) = MIN( ztsmlt(ji) , sist_1d(ji) ) 317 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( sist_1d(ji) - zts_old(ji) ) 318 qla_ice_1d(ji) = qla_ice_1d(ji) + dqla_ice_1d(ji) * ( sist_1d(ji) - zts_old(ji) ) 319 zfcsu(ji) = zksndh(ji) * ( ztbif(ji) - sist_1d(ji) ) 320 END DO 321 ELSE 322 DO ji = kideb, kiut 323 sist_1d(ji) = MIN( ztsmlt(ji) , sist_1d(ji) ) 324 zfcsu(ji) = zksndh(ji) * ( ztbif(ji) - sist_1d(ji) ) 325 END DO 326 ENDIF 327 328 ! 5.2. Calculate available heat for surface ablation. 329 !--------------------------------------------------------------------- 330 331 DO ji = kideb, kiut 332 zfnet(ji) = qns_ice_1d(ji) + z1mi0(ji) * qsr_ice_1d(ji) + zfcsu(ji) 333 zfnet(ji) = MAX( zzero , zfnet(ji) ) 334 zfnet(ji) = zfnet(ji) * MAX( zzero , SIGN( zone , sist_1d(ji) - ztsmlt(ji) ) ) 335 END DO 336 337 !------------------------------------------------------------------------- 338 ! 6. Calculate changes in internal temperature due to vertical diffusion 339 ! processes. The evolution of this temperature is governed by the one- 340 ! dimensionnal heat-diffusion equation. 341 ! Given the temperature tbif(1/2/3), at time m we solve a set 342 ! of finite difference equations to obtain new tempe. Each tempe is coupled 343 ! to the temp. immediatly above and below by heat conduction terms. Thus 344 ! we have a set of equations of the form A * T = B, where A is a tridiagonal 345 ! matrix, T a vector whose components are the unknown new temp. 346 !------------------------------------------------------------------------- 354 347 355 356 357 358 359 360 361 362 363 364 348 !--parameter for the numerical methode use to solve the heat-diffusion equation 349 !- implicit, explicit or Crank-Nicholson 350 zumsb = 1.0 - sbeta 351 DO ji = kideb, kiut 352 zidsn(ji) = MAX ( zzero, SIGN( zone, hsndif - h_snow_1d(ji) ) ) 353 z1midsn(ji) = 1.0 - zidsn(ji) 354 zihic = MAX ( zzero, SIGN( zone, hicdif - h_ice_1d(ji) ) ) 355 zidsnic(ji) = zidsn(ji) * zihic 356 zfcsudt(ji) = zfcsu(ji) * rdt_ice 357 END DO 365 358 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 359 DO ji = kideb, kiut 360 361 ! 6.1 Calculate intermediate variables. 362 !---------------------------------------- 363 364 !--conductivity at the snow surface 365 zkhsn = 2.0 * zksn(ji) * rdt_ice / rcpsn 366 !--conductivity at the ice surface 367 zkhic = 4.0 * zkic(ji) * rdt_ice / MAX( h_ice_1d(ji) * h_ice_1d(ji) * rcpic , epsi20 ) 368 !--conductivity at the snow/ice interface 369 zkint = 4.0 * zksn(ji) * zkic(ji) & 370 & / ( zksn(ji) * h_ice_1d(ji) + 2.0 * zkic(ji) * h_snow_1d(ji) * z1midsn(ji)) 371 zkhsnint = zkint * rdt_ice / rcpsn 372 zkhicint = zkint * 2.0 * rdt_ice / MAX( h_ice_1d(ji) * rcpic , epsi20 ) 373 374 ! 6.2. Fulfill the linear system matrix. 375 !----------------------------------------- 383 376 !$$$ zplediag(1) = 1 + sbeta * z1midsn(ji) * ( zkhsn + zkhsnint ) 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 377 zplediag(1) = zidsn(ji) + z1midsn(ji) * h_snow_1d(ji) & 378 & + sbeta * z1midsn(ji) * zkhsnint 379 zplediag(2) = 1 + sbeta * ( z1midsn(ji) * zkhicint + zkhic ) 380 zplediag(3) = 1 + 3.0 * sbeta * zkhic 381 382 zsubdiag(1) = 0.e0 383 zsubdiag(2) = -1.e0 * z1midsn(ji) * sbeta * zkhicint 384 zsubdiag(3) = -1.e0 * sbeta * zkhic 385 386 zsupdiag(1) = -1.e0 * z1midsn(ji) * sbeta * zkhsnint 387 zsupdiag(2) = zsubdiag(3) 388 zsupdiag(3) = 0.e0 389 390 ! 6.3. Fulfill the idependent term vector. 391 !------------------------------------------- 399 392 400 393 !$$$ zsmbr(1) = zidsn(ji) * sist_1d(ji) + z1midsn(ji) * & … … 402 395 !$$$ & - zumsb * ( zkhsn * tbif_1d(ji,1) 403 396 !$$$ & + zkhsnint * ( tbif_1d(ji,1) - tbif_1d(ji,2) ) ) ) 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 397 zsmbr(1) = zidsn(ji) * sist_1d(ji) + z1midsn(ji) * & 398 & ( h_snow_1d(ji) * tbif_1d(ji,1) - ( zfcsudt(ji) / rcpsn ) & 399 & - zumsb * zkhsnint * ( tbif_1d(ji,1) - tbif_1d(ji,2) ) ) 400 401 zsmbr(2) = tbif_1d(ji,2) & 402 & - zidsn(ji) * ( 1.0 - zidsnic(ji) ) & 403 & * ( zfcsudt(ji) / MAX( h_ice_1d(ji) * rcpic , epsi20 ) ) & 404 & + zumsb * ( zkhicint * ( tbif_1d(ji,1) - tbif_1d(ji,2) ) & 405 & - zkhic * ( tbif_1d(ji,2) - tbif_1d(ji,3) ) ) 406 407 zsmbr(3) = tbif_1d(ji,3) & 408 & + zkhic * ( 2.0 * tfu_1d(ji) & 409 & + zumsb * ( tbif_1d(ji,2) - 3.0 * tbif_1d(ji,3) ) ) 410 411 ! 6.4. Solve the system (Gauss elimination method). 412 !---------------------------------------------------- 413 414 zpiv1 = zsubdiag(2) / zplediag(1) 415 zb2 = zplediag(2) - zpiv1 * zsupdiag(1) 416 zd2 = zsmbr(2) - zpiv1 * zsmbr(1) 417 418 zpiv2 = zsubdiag(3) / zb2 419 zb3 = zplediag(3) - zpiv2 * zsupdiag(2) 420 zd3 = zsmbr(3) - zpiv2 * zd2 421 422 tbif_1d(ji,3) = zd3 / zb3 423 tbif_1d(ji,2) = ( zd2 - zsupdiag(2) * tbif_1d(ji,3) ) / zb2 424 tbif_1d(ji,1) = ( zsmbr(1) - zsupdiag(1) * tbif_1d(ji,2) ) / zplediag(1) 425 426 !--- taking into account the particular case of zidsnic(ji) = 1 427 ztint = ( zkic(ji) * h_snow_1d(ji) * tfu_1d (ji) & 428 & + zksn(ji) * h_ice_1d(ji) * sist_1d(ji) ) & 429 & / ( zkic(ji) * h_snow_1d(ji) + zksn(ji) * h_ice_1d(ji) ) 430 431 tbif_1d(ji,1) = ( 1.0 - zidsnic(ji) ) * tbif_1d(ji,1) & 432 & + zidsnic(ji) * ( ztint + sist_1d(ji) ) / 2.0 433 tbif_1d(ji,2) = ( 1.0 - zidsnic(ji) ) * tbif_1d(ji,2) & 434 & + zidsnic(ji) * ( 3.0 * ztint + tfu_1d(ji) ) / 4.0 435 tbif_1d(ji,3) = ( 1.0 - zidsnic(ji) ) * tbif_1d(ji,3) & 436 & + zidsnic(ji) * ( ztint + 3.0 * tfu_1d(ji) ) / 4.0 437 END DO 445 438 446 447 448 439 !---------------------------------------------------------------------- 440 ! 9. Take into account surface ablation and bottom accretion-ablation.| 441 !---------------------------------------------------------------------- 449 442 450 451 452 453 454 455 456 457 443 !---Snow accumulation in one thermodynamic time step 444 zsprecip(kideb:kiut) = sprecip_1d(kideb:kiut) * rdt_ice / rhosn 445 446 447 DO ji = kideb, kiut 448 449 ! 9.1. Surface ablation and update of snow thickness and qstbif_1d 450 !-------------------------------------------------------------------- 458 451 459 452 !-------------------------------------------------------------------------- … … 759 752 qstbif_1d(ji) = zdrfrl2 * qstbif_1d(ji) 760 753 frld_1d(ji) = zfrld_1d(ji) 761 762 763 764 754 ! 755 END DO 756 ! 757 END SUBROUTINE lim_thd_zdf_2 765 758 766 759 #else … … 768 761 !! Default Option NO sea-ice model 769 762 !!---------------------------------------------------------------------- 770 CONTAINS771 SUBROUTINE lim_thd_zdf_2 ! Empty routine772 END SUBROUTINE lim_thd_zdf_2773 763 #endif 774 764
Note: See TracChangeset
for help on using the changeset viewer.
