- Timestamp:
- 2015-07-10T13:28:53+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4765_CNRS_agrif/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r4765 r5581 19 19 USE phycst ! physical constants (ocean directory) 20 20 USE ice ! LIM-3 variables 21 USE par_ice ! LIM-3 parameters22 21 USE thd_ice ! LIM-3: thermodynamics 23 22 USE in_out_manager ! I/O manager … … 25 24 USE wrk_nemo ! work arrays 26 25 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 27 USE cpl_oasis3, ONLY : lk_cpl28 26 29 27 IMPLICIT NONE … … 32 30 PUBLIC lim_thd_dif ! called by lim_thd 33 31 34 REAL(wp) :: epsi10 = 1.e-10_wp !35 32 !!---------------------------------------------------------------------- 36 33 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) … … 75 72 !! 76 73 !! ** Inputs / Ouputs : (global commons) 77 !! surface temperature : t_su_ b78 !! ice/snow temperatures : t_i_ b, t_s_b79 !! ice salinities : s_i_ b74 !! surface temperature : t_su_1d 75 !! ice/snow temperatures : t_i_1d, t_s_1d 76 !! ice salinities : s_i_1d 80 77 !! number of layers in the ice/snow: nlay_i, nlay_s 81 78 !! profile of the ice/snow layers : z_i, z_s 82 !! total ice/snow thickness : ht_i_ b, ht_s_b79 !! total ice/snow thickness : ht_i_1d, ht_s_1d 83 80 !! 84 81 !! ** External : … … 98 95 INTEGER :: ii, ij ! temporary dummy loop index 99 96 INTEGER :: numeq ! current reference number of equation 100 INTEGER :: layer! vertical dummy loop index97 INTEGER :: jk ! vertical dummy loop index 101 98 INTEGER :: nconv ! number of iterations in iterative procedure 102 99 INTEGER :: minnumeqmin, maxnumeqmax 100 103 101 INTEGER, POINTER, DIMENSION(:) :: numeqmin ! reference number of top equation 104 102 INTEGER, POINTER, DIMENSION(:) :: numeqmax ! reference number of bottom equation 105 INTEGER, POINTER, DIMENSION(:) :: isnow ! switch for presence (1) or absence (0) of snow103 106 104 REAL(wp) :: zg1s = 2._wp ! for the tridiagonal system 107 105 REAL(wp) :: zg1 = 2._wp ! 108 106 REAL(wp) :: zgamma = 18009._wp ! for specific heat 109 107 REAL(wp) :: zbeta = 0.117_wp ! for thermal conductivity (could be 0.13) 110 REAL(wp) :: zraext_s = 1 .e+8_wp! extinction coefficient of radiation in the snow108 REAL(wp) :: zraext_s = 10._wp ! extinction coefficient of radiation in the snow 111 109 REAL(wp) :: zkimin = 0.10_wp ! minimum ice thermal conductivity 112 110 REAL(wp) :: ztsu_err = 1.e-5_wp ! range around which t_su is considered as 0°C 113 111 REAL(wp) :: ztmelt_i ! ice melting temperature 114 112 REAL(wp) :: zerritmax ! current maximal error on temperature 115 REAL(wp), POINTER, DIMENSION(:) :: ztfs ! ice melting point 116 REAL(wp), POINTER, DIMENSION(:) :: ztsuold ! old surface temperature (before the iterative procedure ) 117 REAL(wp), POINTER, DIMENSION(:) :: ztsuoldit ! surface temperature at previous iteration 118 REAL(wp), POINTER, DIMENSION(:) :: zh_i ! ice layer thickness 119 REAL(wp), POINTER, DIMENSION(:) :: zh_s ! snow layer thickness 120 REAL(wp), POINTER, DIMENSION(:) :: zfsw ! solar radiation absorbed at the surface 121 REAL(wp), POINTER, DIMENSION(:) :: zf ! surface flux function 122 REAL(wp), POINTER, DIMENSION(:) :: dzf ! derivative of the surface flux function 123 REAL(wp), POINTER, DIMENSION(:) :: zerrit ! current error on temperature 124 REAL(wp), POINTER, DIMENSION(:) :: zdifcase ! case of the equation resolution (1->4) 125 REAL(wp), POINTER, DIMENSION(:) :: zftrice ! solar radiation transmitted through the ice 126 REAL(wp), POINTER, DIMENSION(:) :: zihic, zhsu 127 REAL(wp), POINTER, DIMENSION(:,:) :: ztcond_i ! Ice thermal conductivity 128 REAL(wp), POINTER, DIMENSION(:,:) :: zradtr_i ! Radiation transmitted through the ice 129 REAL(wp), POINTER, DIMENSION(:,:) :: zradab_i ! Radiation absorbed in the ice 130 REAL(wp), POINTER, DIMENSION(:,:) :: zkappa_i ! Kappa factor in the ice 131 REAL(wp), POINTER, DIMENSION(:,:) :: ztiold ! Old temperature in the ice 132 REAL(wp), POINTER, DIMENSION(:,:) :: zeta_i ! Eta factor in the ice 133 REAL(wp), POINTER, DIMENSION(:,:) :: ztitemp ! Temporary temperature in the ice to check the convergence 134 REAL(wp), POINTER, DIMENSION(:,:) :: zspeche_i ! Ice specific heat 135 REAL(wp), POINTER, DIMENSION(:,:) :: z_i ! Vertical cotes of the layers in the ice 136 REAL(wp), POINTER, DIMENSION(:,:) :: zradtr_s ! Radiation transmited through the snow 137 REAL(wp), POINTER, DIMENSION(:,:) :: zradab_s ! Radiation absorbed in the snow 138 REAL(wp), POINTER, DIMENSION(:,:) :: zkappa_s ! Kappa factor in the snow 139 REAL(wp), POINTER, DIMENSION(:,:) :: zeta_s ! Eta factor in the snow 140 REAL(wp), POINTER, DIMENSION(:,:) :: ztstemp ! Temporary temperature in the snow to check the convergence 141 REAL(wp), POINTER, DIMENSION(:,:) :: ztsold ! Temporary temperature in the snow 142 REAL(wp), POINTER, DIMENSION(:,:) :: z_s ! Vertical cotes of the layers in the snow 143 REAL(wp), POINTER, DIMENSION(:,:) :: zindterm ! Independent term 144 REAL(wp), POINTER, DIMENSION(:,:) :: zindtbis ! temporary independent term 145 REAL(wp), POINTER, DIMENSION(:,:) :: zdiagbis 146 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrid ! tridiagonal system terms 113 REAL(wp) :: zhsu 114 115 REAL(wp), POINTER, DIMENSION(:) :: isnow ! switch for presence (1) or absence (0) of snow 116 REAL(wp), POINTER, DIMENSION(:) :: ztsub ! old surface temperature (before the iterative procedure ) 117 REAL(wp), POINTER, DIMENSION(:) :: ztsubit ! surface temperature at previous iteration 118 REAL(wp), POINTER, DIMENSION(:) :: zh_i ! ice layer thickness 119 REAL(wp), POINTER, DIMENSION(:) :: zh_s ! snow layer thickness 120 REAL(wp), POINTER, DIMENSION(:) :: zfsw ! solar radiation absorbed at the surface 121 REAL(wp), POINTER, DIMENSION(:) :: zqns_ice_b ! solar radiation absorbed at the surface 122 REAL(wp), POINTER, DIMENSION(:) :: zf ! surface flux function 123 REAL(wp), POINTER, DIMENSION(:) :: dzf ! derivative of the surface flux function 124 REAL(wp), POINTER, DIMENSION(:) :: zerrit ! current error on temperature 125 REAL(wp), POINTER, DIMENSION(:) :: zdifcase ! case of the equation resolution (1->4) 126 REAL(wp), POINTER, DIMENSION(:) :: zftrice ! solar radiation transmitted through the ice 127 REAL(wp), POINTER, DIMENSION(:) :: zihic 128 129 REAL(wp), POINTER, DIMENSION(:,:) :: ztcond_i ! Ice thermal conductivity 130 REAL(wp), POINTER, DIMENSION(:,:) :: zradtr_i ! Radiation transmitted through the ice 131 REAL(wp), POINTER, DIMENSION(:,:) :: zradab_i ! Radiation absorbed in the ice 132 REAL(wp), POINTER, DIMENSION(:,:) :: zkappa_i ! Kappa factor in the ice 133 REAL(wp), POINTER, DIMENSION(:,:) :: ztib ! Old temperature in the ice 134 REAL(wp), POINTER, DIMENSION(:,:) :: zeta_i ! Eta factor in the ice 135 REAL(wp), POINTER, DIMENSION(:,:) :: ztitemp ! Temporary temperature in the ice to check the convergence 136 REAL(wp), POINTER, DIMENSION(:,:) :: zspeche_i ! Ice specific heat 137 REAL(wp), POINTER, DIMENSION(:,:) :: z_i ! Vertical cotes of the layers in the ice 138 REAL(wp), POINTER, DIMENSION(:,:) :: zradtr_s ! Radiation transmited through the snow 139 REAL(wp), POINTER, DIMENSION(:,:) :: zradab_s ! Radiation absorbed in the snow 140 REAL(wp), POINTER, DIMENSION(:,:) :: zkappa_s ! Kappa factor in the snow 141 REAL(wp), POINTER, DIMENSION(:,:) :: zeta_s ! Eta factor in the snow 142 REAL(wp), POINTER, DIMENSION(:,:) :: ztstemp ! Temporary temperature in the snow to check the convergence 143 REAL(wp), POINTER, DIMENSION(:,:) :: ztsb ! Temporary temperature in the snow 144 REAL(wp), POINTER, DIMENSION(:,:) :: z_s ! Vertical cotes of the layers in the snow 145 REAL(wp), POINTER, DIMENSION(:,:) :: zindterm ! 'Ind'ependent term 146 REAL(wp), POINTER, DIMENSION(:,:) :: zindtbis ! Temporary 'ind'ependent term 147 REAL(wp), POINTER, DIMENSION(:,:) :: zdiagbis ! Temporary 'dia'gonal term 148 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrid ! Tridiagonal system terms 149 147 150 ! diag errors on heat 148 REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini 149 REAL(wp) :: zhfx_err 151 REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini, zhfx_err 152 153 ! Mono-category 154 REAL(wp) :: zepsilon ! determines thres. above which computation of G(h) is done 155 REAL(wp) :: zratio_s ! dummy factor 156 REAL(wp) :: zratio_i ! dummy factor 157 REAL(wp) :: zh_thres ! thickness thres. for G(h) computation 158 REAL(wp) :: zhe ! dummy factor 159 REAL(wp) :: zkimean ! mean sea ice thermal conductivity 160 REAL(wp) :: zfac ! dummy factor 161 REAL(wp) :: zihe ! dummy factor 162 REAL(wp) :: zheshth ! dummy factor 163 164 REAL(wp), POINTER, DIMENSION(:) :: zghe ! G(he), th. conduct enhancement factor, mono-cat 165 150 166 !!------------------------------------------------------------------ 151 167 ! 152 CALL wrk_alloc( jpij, numeqmin, numeqmax , isnow)153 CALL wrk_alloc( jpij, ztfs, ztsuold, ztsuoldit, zh_i, zh_s, zfsw )154 CALL wrk_alloc( jpij, zf, dzf, z errit, zdifcase, zftrice, zihic, zhsu)155 CALL wrk_alloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztiold, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0)156 CALL wrk_alloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsold, zeta_s, ztstemp, z_s, kjstart=0)157 CALL wrk_alloc( jpij, jkmax+2, zindterm, zindtbis, zdiagbis )158 CALL wrk_alloc( jpij, jkmax+2,3, ztrid )159 160 CALL wrk_alloc( jpij, zdq, zq_ini )168 CALL wrk_alloc( jpij, numeqmin, numeqmax ) 169 CALL wrk_alloc( jpij, isnow, ztsub, ztsubit, zh_i, zh_s, zfsw ) 170 CALL wrk_alloc( jpij, zf, dzf, zqns_ice_b, zerrit, zdifcase, zftrice, zihic, zghe ) 171 CALL wrk_alloc( jpij,nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0 ) 172 CALL wrk_alloc( jpij,nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart=0 ) 173 CALL wrk_alloc( jpij,nlay_i+3, zindterm, zindtbis, zdiagbis ) 174 CALL wrk_alloc( jpij,nlay_i+3,3, ztrid ) 175 176 CALL wrk_alloc( jpij, zdq, zq_ini, zhfx_err ) 161 177 162 178 ! --- diag error on heat diffusion - PART 1 --- ! 163 179 zdq(:) = 0._wp ; zq_ini(:) = 0._wp 164 180 DO ji = kideb, kiut 165 zq_ini(ji) = ( SUM( q_i_ b(ji,1:nlay_i) ) * ht_i_b(ji) / REAL( nlay_i )+ &166 & SUM( q_s_ b(ji,1:nlay_s) ) * ht_s_b(ji) / REAL( nlay_s ))181 zq_ini(ji) = ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i + & 182 & SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s ) 167 183 END DO 168 184 … … 170 186 ! 1) Initialization ! 171 187 !------------------------------------------------------------------------------! 172 ! clem clean: replace just ztfs by rtt173 188 DO ji = kideb , kiut 174 ! is there snow or not 175 isnow(ji)= NINT( 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_b(ji) ) ) ) 176 ! surface temperature of fusion 177 ztfs(ji) = REAL( isnow(ji) ) * rtt + REAL( 1 - isnow(ji) ) * rtt 189 isnow(ji)= 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) ) ! is there snow or not 178 190 ! layer thickness 179 zh_i(ji) = ht_i_ b(ji) / REAL( nlay_i )180 zh_s(ji) = ht_s_ b(ji) / REAL( nlay_s )191 zh_i(ji) = ht_i_1d(ji) * r1_nlay_i 192 zh_s(ji) = ht_s_1d(ji) * r1_nlay_s 181 193 END DO 182 194 … … 188 200 z_i(:,0) = 0._wp ! vert. coord. of the up. lim. of the 1st ice layer 189 201 190 DO layer= 1, nlay_s ! vert. coord of the up. lim. of the layer-th snow layer191 DO ji = kideb , kiut 192 z_s(ji, layer) = z_s(ji,layer-1) + ht_s_b(ji) / REAL( nlay_s )193 END DO 194 END DO 195 196 DO layer= 1, nlay_i ! vert. coord of the up. lim. of the layer-th ice layer197 DO ji = kideb , kiut 198 z_i(ji, layer) = z_i(ji,layer-1) + ht_i_b(ji) / REAL( nlay_i )202 DO jk = 1, nlay_s ! vert. coord of the up. lim. of the layer-th snow layer 203 DO ji = kideb , kiut 204 z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) * r1_nlay_s 205 END DO 206 END DO 207 208 DO jk = 1, nlay_i ! vert. coord of the up. lim. of the layer-th ice layer 209 DO ji = kideb , kiut 210 z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) * r1_nlay_i 199 211 END DO 200 212 END DO 201 213 ! 202 214 !------------------------------------------------------------------------------| 203 ! 2) Radiation s|215 ! 2) Radiation | 204 216 !------------------------------------------------------------------------------| 205 217 ! … … 214 226 ! zftrice = io.qsr_ice is below the surface 215 227 ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice 216 228 ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover 229 zhsu = 0.1_wp ! threshold for the computation of i0 217 230 DO ji = kideb , kiut 218 231 ! switches 219 isnow(ji) = NINT( 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_b(ji) ) ))232 isnow(ji) = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) ) 220 233 ! hs > 0, isnow = 1 221 zhsu (ji) = hnzst ! threshold for the computation of i0 222 zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_b(ji) / zhsu(ji) ) ) 223 224 i0(ji) = REAL( 1 - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) ) 225 !fr1_i0_1d = i0 for a thin ice surface 226 !fr1_i0_2d = i0 for a thick ice surface 227 ! a function of the cloud cover 228 ! 229 !i0(ji) = (1.0-FLOAT(isnow(ji)))*3.0/(100*ht_s_b(ji)+10.0) 230 !formula used in Cice 234 zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) / zhsu ) ) 235 236 i0(ji) = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) ) 231 237 END DO 232 238 … … 236 242 !------------------------------------------------------- 237 243 DO ji = kideb , kiut 238 zfsw (ji) = qsr_ice_1d(ji) * ( 1 - i0(ji) ) ! Shortwave radiation absorbed at surface 239 zftrice(ji) = qsr_ice_1d(ji) * i0(ji) ! Solar radiation transmitted below the surface layer 240 dzf (ji) = dqns_ice_1d(ji) ! derivative of incoming nonsolar flux 244 zfsw (ji) = qsr_ice_1d(ji) * ( 1 - i0(ji) ) ! Shortwave radiation absorbed at surface 245 zftrice(ji) = qsr_ice_1d(ji) * i0(ji) ! Solar radiation transmitted below the surface layer 246 dzf (ji) = dqns_ice_1d(ji) ! derivative of incoming nonsolar flux 247 zqns_ice_b(ji) = qns_ice_1d(ji) ! store previous qns_ice_1d value 241 248 END DO 242 249 … … 249 256 END DO 250 257 251 DO layer= 1, nlay_s ! Radiation through snow258 DO jk = 1, nlay_s ! Radiation through snow 252 259 DO ji = kideb, kiut 253 260 ! ! radiation transmitted below the layer-th snow layer 254 zradtr_s(ji, layer) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,layer) ) ) )261 zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,jk) ) ) ) 255 262 ! ! radiation absorbed by the layer-th snow layer 256 zradab_s(ji, layer) = zradtr_s(ji,layer-1) - zradtr_s(ji,layer)263 zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk) 257 264 END DO 258 265 END DO 259 266 260 267 DO ji = kideb, kiut ! ice initialization 261 zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * REAL( isnow(ji) ) + zftrice(ji) * REAL( 1- isnow(ji) )262 END DO 263 264 DO layer= 1, nlay_i ! Radiation through ice268 zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * isnow(ji) + zftrice(ji) * ( 1._wp - isnow(ji) ) 269 END DO 270 271 DO jk = 1, nlay_i ! Radiation through ice 265 272 DO ji = kideb, kiut 266 273 ! ! radiation transmitted below the layer-th ice layer 267 zradtr_i(ji, layer) = zradtr_i(ji,0) * EXP( - kappa_i * ( MAX ( 0._wp , z_i(ji,layer) ) ) )274 zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * ( MAX ( 0._wp , z_i(ji,jk) ) ) ) 268 275 ! ! radiation absorbed by the layer-th ice layer 269 zradab_i(ji, layer) = zradtr_i(ji,layer-1) - zradtr_i(ji,layer)276 zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk) 270 277 END DO 271 278 END DO 272 279 273 280 DO ji = kideb, kiut ! Radiation transmitted below the ice 274 !!!ftr_ice_1d(ji) = ftr_ice_1d(ji) + iatte_1d(ji) * zradtr_i(ji,nlay_i) * a_i_b(ji) / at_i_b(ji) ! clem modif275 281 ftr_ice_1d(ji) = zradtr_i(ji,nlay_i) 276 282 END DO 277 283 278 !279 284 !------------------------------------------------------------------------------| 280 285 ! 3) Iterative procedure begins | … … 282 287 ! 283 288 DO ji = kideb, kiut ! Old surface temperature 284 ztsu old (ji) = t_su_b(ji) ! temperature at the beg of iter pr.285 ztsu oldit(ji) = t_su_b(ji) ! temperature at the previous iter286 t_su_ b (ji) = MIN( t_su_b(ji), ztfs(ji) - ztsu_err )! necessary287 zerrit (ji) = 1000._wp! initial value of error288 END DO 289 290 DO layer= 1, nlay_s ! Old snow temperature291 DO ji = kideb , kiut 292 zts old(ji,layer) = t_s_b(ji,layer)293 END DO 294 END DO 295 296 DO layer= 1, nlay_i ! Old ice temperature297 DO ji = kideb , kiut 298 zti old(ji,layer) = t_i_b(ji,layer)289 ztsub (ji) = t_su_1d(ji) ! temperature at the beg of iter pr. 290 ztsubit(ji) = t_su_1d(ji) ! temperature at the previous iter 291 t_su_1d(ji) = MIN( t_su_1d(ji), rt0 - ztsu_err ) ! necessary 292 zerrit (ji) = 1000._wp ! initial value of error 293 END DO 294 295 DO jk = 1, nlay_s ! Old snow temperature 296 DO ji = kideb , kiut 297 ztsb(ji,jk) = t_s_1d(ji,jk) 298 END DO 299 END DO 300 301 DO jk = 1, nlay_i ! Old ice temperature 302 DO ji = kideb , kiut 303 ztib(ji,jk) = t_i_1d(ji,jk) 299 304 END DO 300 305 END DO … … 303 308 zerritmax = 1000._wp ! maximal value of error on all points 304 309 305 DO WHILE ( zerritmax > maxer_i_thd .AND. nconv < nconv_i_thd)310 DO WHILE ( zerritmax > rn_terr_dif .AND. nconv < nn_conv_dif ) 306 311 ! 307 312 nconv = nconv + 1 … … 311 316 !------------------------------------------------------------------------------| 312 317 ! 313 IF( thcon_i_swi== 0 ) THEN ! Untersteiner (1964) formula314 DO ji = kideb , kiut 315 ztcond_i(ji,0) = rcdic + zbeta*s_i_b(ji,1) / MIN(-epsi10,t_i_b(ji,1)-rtt)316 ztcond_i(ji,0) = MAX(ztcond_i(ji,0),zkimin)317 END DO 318 DO layer= 1, nlay_i-1318 IF( nn_ice_thcon == 0 ) THEN ! Untersteiner (1964) formula 319 DO ji = kideb , kiut 320 ztcond_i(ji,0) = rcdic + zbeta * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) 321 ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin ) 322 END DO 323 DO jk = 1, nlay_i-1 319 324 DO ji = kideb , kiut 320 ztcond_i(ji, layer) = rcdic + zbeta*( s_i_b(ji,layer) + s_i_b(ji,layer+1) ) / &321 MIN(-2.0_wp * epsi10, t_i_ b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt)322 ztcond_i(ji, layer) = MAX(ztcond_i(ji,layer),zkimin)325 ztcond_i(ji,jk) = rcdic + zbeta * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) / & 326 MIN(-2.0_wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0) 327 ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin ) 323 328 END DO 324 329 END DO 325 330 ENDIF 326 331 327 IF( thcon_i_swi== 1 ) THEN ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T328 DO ji = kideb , kiut 329 ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_ b(ji,1) / MIN( -epsi10, t_i_b(ji,1)-rtt) &330 & - 0.011_wp * ( t_i_ b(ji,1) - rtt)332 IF( nn_ice_thcon == 1 ) THEN ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T 333 DO ji = kideb , kiut 334 ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) & 335 & - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) 331 336 ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin ) 332 337 END DO 333 DO layer= 1, nlay_i-1338 DO jk = 1, nlay_i-1 334 339 DO ji = kideb , kiut 335 ztcond_i(ji, layer) = rcdic +&336 & 0.09 0_wp * ( s_i_b(ji,layer) + s_i_b(ji,layer+1) )&337 & / MIN( -2.0_wp * epsi10, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt) &338 & - 0.0055_wp * ( t_i_b(ji,layer) + t_i_b(ji,layer+1) - 2.0*rtt)339 ztcond_i(ji, layer) = MAX( ztcond_i(ji,layer), zkimin )340 ztcond_i(ji,jk) = rcdic + & 341 & 0.09_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) & 342 & / MIN( -2._wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0 ) & 343 & - 0.0055_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0 * rt0 ) 344 ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin ) 340 345 END DO 341 346 END DO 342 347 DO ji = kideb , kiut 343 ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_ b(ji,nlay_i) / MIN(-epsi10,t_bo_b(ji)-rtt) &344 & - 0.011_wp * ( t_bo_ b(ji) - rtt)348 ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) & 349 & - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 345 350 ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin ) 346 351 END DO 347 352 ENDIF 348 ! 349 !------------------------------------------------------------------------------| 350 ! 5) kappa factors | 351 !------------------------------------------------------------------------------| 352 ! 353 354 ! 355 !------------------------------------------------------------------------------| 356 ! 5) G(he) - enhancement of thermal conductivity in mono-category case | 357 !------------------------------------------------------------------------------| 358 ! 359 ! Computation of effective thermal conductivity G(h) 360 ! Used in mono-category case only to simulate an ITD implicitly 361 ! Fichefet and Morales Maqueda, JGR 1997 362 363 zghe(:) = 1._wp 364 365 SELECT CASE ( nn_monocat ) 366 367 CASE (1,3) ! LIM3 368 369 zepsilon = 0.1_wp 370 zh_thres = EXP( 1._wp ) * zepsilon * 0.5_wp 371 372 DO ji = kideb, kiut 373 374 ! Mean sea ice thermal conductivity 375 zkimean = SUM( ztcond_i(ji,0:nlay_i) ) / REAL( nlay_i+1, wp ) 376 377 ! Effective thickness he (zhe) 378 zfac = 1._wp / ( rcdsn + zkimean ) 379 zratio_s = rcdsn * zfac 380 zratio_i = zkimean * zfac 381 zhe = zratio_s * ht_i_1d(ji) + zratio_i * ht_s_1d(ji) 382 383 ! G(he) 384 rswitch = MAX( 0._wp , SIGN( 1._wp , zhe - zh_thres ) ) ! =0 if zhe < zh_thres, if > 385 zghe(ji) = ( 1._wp - rswitch ) + rswitch * 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) 386 387 ! Impose G(he) < 2. 388 zghe(ji) = MIN( zghe(ji), 2._wp ) 389 390 END DO 391 392 END SELECT 393 394 ! 395 !------------------------------------------------------------------------------| 396 ! 6) kappa factors | 397 !------------------------------------------------------------------------------| 398 ! 399 !--- Snow 353 400 DO ji = kideb, kiut 354 355 !-- Snow kappa factors 356 zkappa_s(ji,0) = rcdsn / MAX(epsi10,zh_s(ji)) 357 zkappa_s(ji,nlay_s) = rcdsn / MAX(epsi10,zh_s(ji)) 358 END DO 359 360 DO layer = 1, nlay_s-1 361 DO ji = kideb , kiut 362 zkappa_s(ji,layer) = 2.0 * rcdsn / & 363 MAX(epsi10,2.0*zh_s(ji)) 364 END DO 365 END DO 366 367 DO layer = 1, nlay_i-1 368 DO ji = kideb , kiut 369 !-- Ice kappa factors 370 zkappa_i(ji,layer) = 2.0*ztcond_i(ji,layer)/ & 371 MAX(epsi10,2.0*zh_i(ji)) 372 END DO 373 END DO 374 375 DO ji = kideb , kiut 376 zkappa_i(ji,0) = ztcond_i(ji,0)/MAX(epsi10,zh_i(ji)) 377 zkappa_i(ji,nlay_i) = ztcond_i(ji,nlay_i) / MAX(epsi10,zh_i(ji)) 378 !-- Interface 379 zkappa_s(ji,nlay_s) = 2.0*rcdsn*ztcond_i(ji,0)/MAX(epsi10, & 380 (ztcond_i(ji,0)*zh_s(ji) + rcdsn*zh_i(ji))) 381 zkappa_i(ji,0) = zkappa_s(ji,nlay_s)*REAL( isnow(ji) ) & 382 + zkappa_i(ji,0)*REAL( 1 - isnow(ji) ) 383 END DO 384 ! 385 !------------------------------------------------------------------------------| 386 ! 6) Sea ice specific heat, eta factors | 387 !------------------------------------------------------------------------------| 388 ! 389 DO layer = 1, nlay_i 390 DO ji = kideb , kiut 391 ztitemp(ji,layer) = t_i_b(ji,layer) 392 zspeche_i(ji,layer) = cpic + zgamma*s_i_b(ji,layer)/ & 393 MAX((t_i_b(ji,layer)-rtt)*(ztiold(ji,layer)-rtt),epsi10) 394 zeta_i(ji,layer) = rdt_ice / MAX(rhoic*zspeche_i(ji,layer)*zh_i(ji), & 395 epsi10) 396 END DO 397 END DO 398 399 DO layer = 1, nlay_s 400 DO ji = kideb , kiut 401 ztstemp(ji,layer) = t_s_b(ji,layer) 402 zeta_s(ji,layer) = rdt_ice / MAX(rhosn*cpic*zh_s(ji),epsi10) 403 END DO 404 END DO 405 ! 406 !------------------------------------------------------------------------------| 407 ! 7) surface flux computation | 408 !------------------------------------------------------------------------------| 409 ! 410 DO ji = kideb , kiut 411 ! update of the non solar flux according to the update in T_su 412 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_b(ji) - ztsuoldit(ji) ) 413 401 zfac = 1. / MAX( epsi10 , zh_s(ji) ) 402 zkappa_s(ji,0) = zghe(ji) * rcdsn * zfac 403 zkappa_s(ji,nlay_s) = zghe(ji) * rcdsn * zfac 404 END DO 405 406 DO jk = 1, nlay_s-1 407 DO ji = kideb , kiut 408 zkappa_s(ji,jk) = zghe(ji) * 2.0 * rcdsn / MAX( epsi10, 2.0 * zh_s(ji) ) 409 END DO 410 END DO 411 412 !--- Ice 413 DO jk = 1, nlay_i-1 414 DO ji = kideb , kiut 415 zkappa_i(ji,jk) = zghe(ji) * 2.0 * ztcond_i(ji,jk) / MAX( epsi10 , 2.0 * zh_i(ji) ) 416 END DO 417 END DO 418 419 !--- Snow-ice interface 420 DO ji = kideb , kiut 421 zfac = 1./ MAX( epsi10 , zh_i(ji) ) 422 zkappa_i(ji,0) = zghe(ji) * ztcond_i(ji,0) * zfac 423 zkappa_i(ji,nlay_i) = zghe(ji) * ztcond_i(ji,nlay_i) * zfac 424 zkappa_s(ji,nlay_s) = zghe(ji) * zghe(ji) * 2.0 * rcdsn * ztcond_i(ji,0) / & 425 & MAX( epsi10, ( zghe(ji) * ztcond_i(ji,0) * zh_s(ji) + zghe(ji) * rcdsn * zh_i(ji) ) ) 426 zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 427 END DO 428 429 ! 430 !------------------------------------------------------------------------------| 431 ! 7) Sea ice specific heat, eta factors | 432 !------------------------------------------------------------------------------| 433 ! 434 DO jk = 1, nlay_i 435 DO ji = kideb , kiut 436 ztitemp(ji,jk) = t_i_1d(ji,jk) 437 zspeche_i(ji,jk) = cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztib(ji,jk) - rt0 ), epsi10 ) 438 zeta_i(ji,jk) = rdt_ice / MAX( rhoic * zspeche_i(ji,jk) * zh_i(ji), epsi10 ) 439 END DO 440 END DO 441 442 DO jk = 1, nlay_s 443 DO ji = kideb , kiut 444 ztstemp(ji,jk) = t_s_1d(ji,jk) 445 zeta_s(ji,jk) = rdt_ice / MAX( rhosn * cpic * zh_s(ji), epsi10 ) 446 END DO 447 END DO 448 449 ! 450 !------------------------------------------------------------------------------| 451 ! 8) surface flux computation | 452 !------------------------------------------------------------------------------| 453 ! 454 IF ( ln_it_qnsice ) THEN 455 DO ji = kideb , kiut 456 ! update of the non solar flux according to the update in T_su 457 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsubit(ji) ) 458 END DO 459 ENDIF 460 461 ! Update incoming flux 462 DO ji = kideb , kiut 414 463 ! update incoming flux 415 zf(ji) = zfsw(ji) & ! net absorbed solar radiation 416 + qns_ice_1d(ji) ! non solar total flux 417 ! (LWup, LWdw, SH, LH) 418 END DO 419 420 ! 421 !------------------------------------------------------------------------------| 422 ! 8) tridiagonal system terms | 464 zf(ji) = zfsw(ji) & ! net absorbed solar radiation 465 & + qns_ice_1d(ji) ! non solar total flux (LWup, LWdw, SH, LH) 466 END DO 467 468 ! 469 !------------------------------------------------------------------------------| 470 ! 9) tridiagonal system terms | 423 471 !------------------------------------------------------------------------------| 424 472 ! … … 430 478 !!ice interior terms (top equation has the same form as the others) 431 479 432 DO numeq=1, jkmax+2480 DO numeq=1,nlay_i+3 433 481 DO ji = kideb , kiut 434 482 ztrid(ji,numeq,1) = 0. … … 443 491 DO numeq = nlay_s + 2, nlay_s + nlay_i 444 492 DO ji = kideb , kiut 445 layer = numeq - nlay_s - 1 446 ztrid(ji,numeq,1) = - zeta_i(ji,layer)*zkappa_i(ji,layer-1) 447 ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,layer)*(zkappa_i(ji,layer-1) + & 448 zkappa_i(ji,layer)) 449 ztrid(ji,numeq,3) = - zeta_i(ji,layer)*zkappa_i(ji,layer) 450 zindterm(ji,numeq) = ztiold(ji,layer) + zeta_i(ji,layer)* & 451 zradab_i(ji,layer) 493 jk = numeq - nlay_s - 1 494 ztrid(ji,numeq,1) = - zeta_i(ji,jk) * zkappa_i(ji,jk-1) 495 ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) ) 496 ztrid(ji,numeq,3) = - zeta_i(ji,jk) * zkappa_i(ji,jk) 497 zindterm(ji,numeq) = ztib(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk) 452 498 END DO 453 499 ENDDO … … 457 503 !!ice bottom term 458 504 ztrid(ji,numeq,1) = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1) 459 ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,nlay_i)*( zkappa_i(ji,nlay_i)*zg1 & 460 + zkappa_i(ji,nlay_i-1) ) 505 ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) ) 461 506 ztrid(ji,numeq,3) = 0.0 462 zindterm(ji,numeq) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i)* & 463 ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i)*zg1 & 464 * t_bo_b(ji) ) 507 zindterm(ji,numeq) = ztib(ji,nlay_i) + zeta_i(ji,nlay_i) * & 508 & ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 465 509 ENDDO 466 510 467 511 468 512 DO ji = kideb , kiut 469 IF ( ht_s_ b(ji).gt.0.0 ) THEN513 IF ( ht_s_1d(ji) > 0.0 ) THEN 470 514 ! 471 515 !------------------------------------------------------------------------------| … … 475 519 !!snow interior terms (bottom equation has the same form as the others) 476 520 DO numeq = 3, nlay_s + 1 477 layer = numeq - 1 478 ztrid(ji,numeq,1) = - zeta_s(ji,layer)*zkappa_s(ji,layer-1) 479 ztrid(ji,numeq,2) = 1.0 + zeta_s(ji,layer)*( zkappa_s(ji,layer-1) + & 480 zkappa_s(ji,layer) ) 481 ztrid(ji,numeq,3) = - zeta_s(ji,layer)*zkappa_s(ji,layer) 482 zindterm(ji,numeq) = ztsold(ji,layer) + zeta_s(ji,layer)* & 483 zradab_s(ji,layer) 521 jk = numeq - 1 522 ztrid(ji,numeq,1) = - zeta_s(ji,jk) * zkappa_s(ji,jk-1) 523 ztrid(ji,numeq,2) = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) ) 524 ztrid(ji,numeq,3) = - zeta_s(ji,jk)*zkappa_s(ji,jk) 525 zindterm(ji,numeq) = ztsb(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 484 526 END DO 485 527 … … 487 529 IF ( nlay_i.eq.1 ) THEN 488 530 ztrid(ji,nlay_s+2,3) = 0.0 489 zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1)* & 490 t_bo_b(ji) 531 zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) 491 532 ENDIF 492 533 493 IF ( t_su_ b(ji) .LT. rtt) THEN534 IF ( t_su_1d(ji) < rt0 ) THEN 494 535 495 536 !------------------------------------------------------------------------------| … … 501 542 502 543 !!surface equation 503 ztrid(ji,1,1) = 0.0504 ztrid(ji,1,2) = dzf(ji) - zg1s*zkappa_s(ji,0)505 ztrid(ji,1,3) = zg1s*zkappa_s(ji,0)506 zindterm(ji,1) = dzf(ji) *t_su_b(ji)- zf(ji)544 ztrid(ji,1,1) = 0.0 545 ztrid(ji,1,2) = dzf(ji) - zg1s * zkappa_s(ji,0) 546 ztrid(ji,1,3) = zg1s * zkappa_s(ji,0) 547 zindterm(ji,1) = dzf(ji) * t_su_1d(ji) - zf(ji) 507 548 508 549 !!first layer of snow equation 509 ztrid(ji,2,1) = - zkappa_s(ji,0) *zg1s*zeta_s(ji,1)510 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) *(zkappa_s(ji,1) + zkappa_s(ji,0)*zg1s)550 ztrid(ji,2,1) = - zkappa_s(ji,0) * zg1s * zeta_s(ji,1) 551 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 511 552 ztrid(ji,2,3) = - zeta_s(ji,1)* zkappa_s(ji,1) 512 zindterm(ji,2) = zts old(ji,1) + zeta_s(ji,1)*zradab_s(ji,1)553 zindterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) * zradab_s(ji,1) 513 554 514 555 ELSE … … 524 565 !!first layer of snow equation 525 566 ztrid(ji,2,1) = 0.0 526 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + & 527 zkappa_s(ji,0) * zg1s ) 567 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 528 568 ztrid(ji,2,3) = - zeta_s(ji,1)*zkappa_s(ji,1) 529 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * & 530 ( zradab_s(ji,1) + & 531 zkappa_s(ji,0) * zg1s * t_su_b(ji) ) 569 zindterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) * & 570 & ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 532 571 ENDIF 533 572 ELSE … … 537 576 !------------------------------------------------------------------------------| 538 577 ! 539 IF ( t_su_b(ji) .LT. rtt) THEN578 IF ( t_su_1d(ji) < rt0 ) THEN 540 579 ! 541 580 !------------------------------------------------------------------------------| … … 551 590 ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0)*zg1 552 591 ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0)*zg1 553 zindterm(ji,numeqmin(ji)) = dzf(ji)*t_su_ b(ji) - zf(ji)592 zindterm(ji,numeqmin(ji)) = dzf(ji)*t_su_1d(ji) - zf(ji) 554 593 555 594 !!first layer of ice equation 556 595 ztrid(ji,numeqmin(ji)+1,1) = - zkappa_i(ji,0) * zg1 * zeta_i(ji,1) 557 ztrid(ji,numeqmin(ji)+1,2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) & 558 + zkappa_i(ji,0) * zg1 ) 559 ztrid(ji,numeqmin(ji)+1,3) = - zeta_i(ji,1)*zkappa_i(ji,1) 560 zindterm(ji,numeqmin(ji)+1)= ztiold(ji,1) + zeta_i(ji,1)*zradab_i(ji,1) 596 ztrid(ji,numeqmin(ji)+1,2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 597 ztrid(ji,numeqmin(ji)+1,3) = - zeta_i(ji,1) * zkappa_i(ji,1) 598 zindterm(ji,numeqmin(ji)+1)= ztib(ji,1) + zeta_i(ji,1) * zradab_i(ji,1) 561 599 562 600 !!case of only one layer in the ice (surface & ice equations are altered) 563 601 564 IF ( nlay_i.eq.1) THEN602 IF ( nlay_i == 1 ) THEN 565 603 ztrid(ji,numeqmin(ji),1) = 0.0 566 ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0)*2.0 567 ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0)*2.0 568 ztrid(ji,numeqmin(ji)+1,1) = -zkappa_i(ji,0)*2.0*zeta_i(ji,1) 569 ztrid(ji,numeqmin(ji)+1,2) = 1.0 + zeta_i(ji,1)*(zkappa_i(ji,0)*2.0 + & 570 zkappa_i(ji,1)) 604 ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0) * 2.0 605 ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0) * 2.0 606 ztrid(ji,numeqmin(ji)+1,1) = -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1) 607 ztrid(ji,numeqmin(ji)+1,2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) 571 608 ztrid(ji,numeqmin(ji)+1,3) = 0.0 572 609 573 zindterm(ji,numeqmin(ji)+1) = zti old(ji,1) + zeta_i(ji,1)*&574 ( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_b(ji) )610 zindterm(ji,numeqmin(ji)+1) = ztib(ji,1) + zeta_i(ji,1) * & 611 & ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) 575 612 ENDIF 576 613 … … 588 625 !!first layer of ice equation 589 626 ztrid(ji,numeqmin(ji),1) = 0.0 590 ztrid(ji,numeqmin(ji),2) = 1.0 + zeta_i(ji,1)*(zkappa_i(ji,1) + zkappa_i(ji,0)* & 591 zg1) 627 ztrid(ji,numeqmin(ji),2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 592 628 ztrid(ji,numeqmin(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1) 593 zindterm(ji,numeqmin(ji)) = zti old(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) +&594 zkappa_i(ji,0) * zg1 * t_su_b(ji) )629 zindterm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1) * & 630 & ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 595 631 596 632 !!case of only one layer in the ice (surface & ice equations are altered) 597 IF ( nlay_i.eq.1) THEN633 IF ( nlay_i == 1 ) THEN 598 634 ztrid(ji,numeqmin(ji),1) = 0.0 599 ztrid(ji,numeqmin(ji),2) = 1.0 + zeta_i(ji,1)*(zkappa_i(ji,0)*2.0 + & 600 zkappa_i(ji,1)) 635 ztrid(ji,numeqmin(ji),2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) 601 636 ztrid(ji,numeqmin(ji),3) = 0.0 602 zindterm(ji,numeqmin(ji)) = ztiold(ji,1) + zeta_i(ji,1)* & 603 (zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_b(ji)) & 604 + t_su_b(ji)*zeta_i(ji,1)*zkappa_i(ji,0)*2.0 637 zindterm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) & 638 & + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0 605 639 ENDIF 606 640 … … 612 646 ! 613 647 !------------------------------------------------------------------------------| 614 ! 9) tridiagonal system solving|648 ! 10) tridiagonal system solving | 615 649 !------------------------------------------------------------------------------| 616 650 ! … … 621 655 622 656 maxnumeqmax = 0 623 minnumeqmin = jkmax+4657 minnumeqmin = nlay_i+5 624 658 625 659 DO ji = kideb , kiut … … 630 664 END DO 631 665 632 DO layer = minnumeqmin+1, maxnumeqmax 633 DO ji = kideb , kiut 634 numeq = min(max(numeqmin(ji)+1,layer),numeqmax(ji)) 635 zdiagbis(ji,numeq) = ztrid(ji,numeq,2) - ztrid(ji,numeq,1)* & 636 ztrid(ji,numeq-1,3)/zdiagbis(ji,numeq-1) 637 zindtbis(ji,numeq) = zindterm(ji,numeq) - ztrid(ji,numeq,1)* & 638 zindtbis(ji,numeq-1)/zdiagbis(ji,numeq-1) 666 DO jk = minnumeqmin+1, maxnumeqmax 667 DO ji = kideb , kiut 668 numeq = min(max(numeqmin(ji)+1,jk),numeqmax(ji)) 669 zdiagbis(ji,numeq) = ztrid(ji,numeq,2) - ztrid(ji,numeq,1) * ztrid(ji,numeq-1,3) / zdiagbis(ji,numeq-1) 670 zindtbis(ji,numeq) = zindterm(ji,numeq) - ztrid(ji,numeq,1) * zindtbis(ji,numeq-1) / zdiagbis(ji,numeq-1) 639 671 END DO 640 672 END DO … … 642 674 DO ji = kideb , kiut 643 675 ! ice temperatures 644 t_i_b(ji,nlay_i) = zindtbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji)) 645 END DO 646 647 DO numeq = nlay_i + nlay_s + 1, nlay_s + 2, -1 648 DO ji = kideb , kiut 649 layer = numeq - nlay_s - 1 650 t_i_b(ji,layer) = (zindtbis(ji,numeq) - ztrid(ji,numeq,3)* & 651 t_i_b(ji,layer+1))/zdiagbis(ji,numeq) 676 t_i_1d(ji,nlay_i) = zindtbis(ji,numeqmax(ji)) / zdiagbis(ji,numeqmax(ji)) 677 END DO 678 679 DO numeq = nlay_i + nlay_s, nlay_s + 2, -1 680 DO ji = kideb , kiut 681 jk = numeq - nlay_s - 1 682 t_i_1d(ji,jk) = ( zindtbis(ji,numeq) - ztrid(ji,numeq,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,numeq) 652 683 END DO 653 684 END DO … … 655 686 DO ji = kideb , kiut 656 687 ! snow temperatures 657 IF (ht_s_b(ji).GT.0._wp) & 658 t_s_b(ji,nlay_s) = (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) & 659 * t_i_b(ji,1))/zdiagbis(ji,nlay_s+1) & 660 * MAX(0.0,SIGN(1.0,ht_s_b(ji))) 688 IF (ht_s_1d(ji) > 0._wp) & 689 t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) & 690 & / zdiagbis(ji,nlay_s+1) * MAX( 0.0, SIGN( 1.0, ht_s_1d(ji) ) ) 661 691 662 692 ! surface temperature 663 isnow(ji) = NINT( 1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_b(ji) ) ))664 ztsu oldit(ji) = t_su_b(ji)665 IF( t_su_ b(ji) < ztfs(ji)) &666 t_su_ b(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_b(ji,1)&667 & + REAL( 1 - isnow(ji) )*t_i_b(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))693 isnow(ji) = 1._wp - MAX( 0._wp , SIGN( 1._wp , -ht_s_1d(ji) ) ) 694 ztsubit(ji) = t_su_1d(ji) 695 IF( t_su_1d(ji) < rt0 ) & 696 t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3) * & 697 & ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji)) 668 698 END DO 669 699 ! 670 700 !-------------------------------------------------------------------------- 671 ! 1 0) Has the scheme converged ?, end of the iterative procedure |701 ! 11) Has the scheme converged ?, end of the iterative procedure | 672 702 !-------------------------------------------------------------------------- 673 703 ! 674 704 ! check that nowhere it has started to melt 675 ! zerrit(ji) is a measure of error, it has to be under maxer_i_thd676 DO ji = kideb , kiut 677 t_su_ b(ji) = MAX( MIN( t_su_b(ji) , ztfs(ji)) , 190._wp )678 zerrit(ji) = ABS( t_su_b(ji) - ztsuoldit(ji) )679 END DO 680 681 DO layer= 1, nlay_s682 DO ji = kideb , kiut 683 t_s_ b(ji,layer) = MAX( MIN( t_s_b(ji,layer), rtt), 190._wp )684 zerrit(ji) = MAX(zerrit(ji),ABS(t_s_b(ji,layer) - ztstemp(ji,layer)))685 END DO 686 END DO 687 688 DO layer= 1, nlay_i689 DO ji = kideb , kiut 690 ztmelt_i = -tmut * s_i_b(ji,layer) + rtt691 t_i_ b(ji,layer) = MAX(MIN(t_i_b(ji,layer),ztmelt_i), 190._wp)692 zerrit(ji) = MAX(zerrit(ji),ABS(t_i_b(ji,layer) - ztitemp(ji,layer)))705 ! zerrit(ji) is a measure of error, it has to be under terr_dif 706 DO ji = kideb , kiut 707 t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , 190._wp ) 708 zerrit(ji) = ABS( t_su_1d(ji) - ztsubit(ji) ) 709 END DO 710 711 DO jk = 1, nlay_s 712 DO ji = kideb , kiut 713 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), 190._wp ) 714 zerrit(ji) = MAX( zerrit(ji), ABS( t_s_1d(ji,jk) - ztstemp(ji,jk) ) ) 715 END DO 716 END DO 717 718 DO jk = 1, nlay_i 719 DO ji = kideb , kiut 720 ztmelt_i = -tmut * s_i_1d(ji,jk) + rt0 721 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), 190._wp ) 722 zerrit(ji) = MAX( zerrit(ji), ABS( t_i_1d(ji,jk) - ztitemp(ji,jk) ) ) 693 723 END DO 694 724 END DO … … 704 734 END DO ! End of the do while iterative procedure 705 735 706 IF( ln_ nicep.AND. lwp ) THEN736 IF( ln_icectl .AND. lwp ) THEN 707 737 WRITE(numout,*) ' zerritmax : ', zerritmax 708 738 WRITE(numout,*) ' nconv : ', nconv … … 711 741 ! 712 742 !-------------------------------------------------------------------------! 713 ! 1 1) Fluxes at the interfaces !743 ! 12) Fluxes at the interfaces ! 714 744 !-------------------------------------------------------------------------! 715 745 DO ji = kideb, kiut 716 ! forced mode only : update of latent heat fluxes (sublimation) (always >=0, upward flux)717 IF( .NOT. lk_cpl) qla_ice_1d (ji) = MAX( 0._wp, qla_ice_1d (ji) + dqla_ice_1d(ji) * ( t_su_b(ji) - ztsuold(ji) ) )718 746 ! ! surface ice conduction flux 719 isnow(ji) = NINT( 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_b(ji) ) ))720 fc_su(ji) = - REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * (t_s_b(ji,1) - t_su_b(ji)) &721 & - REAL( 1 - isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_b(ji,1) - t_su_b(ji))747 isnow(ji) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_1d(ji) ) ) 748 fc_su(ji) = - isnow(ji) * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji)) & 749 & - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_1d(ji,1) - t_su_1d(ji)) 722 750 ! ! bottom ice conduction flux 723 fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) ) 724 END DO 751 fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) ) 752 END DO 753 754 ! --- computes sea ice energy of melting compulsory for limthd_dh --- ! 755 CALL lim_thd_enmelt( kideb, kiut ) 756 757 ! --- diagnose the change in non-solar flux due to surface temperature change --- ! 758 IF ( ln_it_qnsice ) THEN 759 DO ji = kideb, kiut 760 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji) - zqns_ice_b(ji) ) * a_i_1d(ji) 761 END DO 762 END IF 763 764 ! --- diag conservation imbalance on heat diffusion - PART 2 --- ! 765 DO ji = kideb, kiut 766 zdq(ji) = - zq_ini(ji) + ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i + & 767 & SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s ) 768 IF( t_su_1d(ji) < rt0 ) THEN ! case T_su < 0degC 769 zhfx_err(ji) = qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice 770 ELSE ! case T_su = 0degC 771 zhfx_err(ji) = fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice 772 ENDIF 773 hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err(ji) * a_i_1d(ji) 774 775 ! total heat that is sent to the ocean (i.e. not used in the heat diffusion equation) 776 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err(ji) * a_i_1d(ji) 777 END DO 725 778 726 779 !----------------------------------------- … … 728 781 !----------------------------------------- 729 782 DO ji = kideb, kiut 730 IF( t_su_ b(ji) < rtt) THEN ! case T_su < 0degC783 IF( t_su_1d(ji) < rt0 ) THEN ! case T_su < 0degC 731 784 hfx_dif_1d(ji) = hfx_dif_1d(ji) + & 732 & ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_ b(ji)733 ELSE ! case T_su = 0degC785 & ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji) 786 ELSE ! case T_su = 0degC 734 787 hfx_dif_1d(ji) = hfx_dif_1d(ji) + & 735 & ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_ b(ji)788 & ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji) 736 789 ENDIF 737 END DO 738 739 ! --- computes sea ice energy of melting compulsory for limthd_dh --- ! 740 CALL lim_thd_enmelt( kideb, kiut ) 741 742 ! --- diag error on heat diffusion - PART 2 --- ! 743 DO ji = kideb, kiut 744 zdq(ji) = - zq_ini(ji) + ( SUM( q_i_b(ji,1:nlay_i) ) * ht_i_b(ji) / REAL( nlay_i ) + & 745 & SUM( q_s_b(ji,1:nlay_s) ) * ht_s_b(ji) / REAL( nlay_s ) ) 746 zhfx_err = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice ) 747 hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err * a_i_b(ji) 748 ! --- correction of qns_ice and surface conduction flux --- ! 749 qns_ice_1d(ji) = qns_ice_1d(ji) - zhfx_err 750 fc_su (ji) = fc_su (ji) - zhfx_err 751 ! --- Heat flux at the ice surface in W.m-2 --- ! 752 ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 753 hfx_in (ii,ij) = hfx_in (ii,ij) + a_i_b(ji) * ( qsr_ice_1d(ji) + qns_ice_1d(ji) ) 754 END DO 755 790 ! correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation 791 hfx_dif_1d(ji) = hfx_dif_1d(ji) - zhfx_err(ji) * a_i_1d(ji) 792 END DO 756 793 ! 757 CALL wrk_dealloc( jpij, numeqmin, numeqmax, isnow ) 758 CALL wrk_dealloc( jpij, ztfs, ztsuold, ztsuoldit, zh_i, zh_s, zfsw ) 759 CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zhsu ) 760 CALL wrk_dealloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, & 761 & ztiold, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 ) 762 CALL wrk_dealloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsold, zeta_s, ztstemp, z_s, kjstart = 0 ) 763 CALL wrk_dealloc( jpij, jkmax+2, zindterm, zindtbis, zdiagbis ) 764 CALL wrk_dealloc( jpij, jkmax+2, 3, ztrid ) 765 CALL wrk_dealloc( jpij, zdq, zq_ini ) 794 CALL wrk_dealloc( jpij, numeqmin, numeqmax ) 795 CALL wrk_dealloc( jpij, isnow, ztsub, ztsubit, zh_i, zh_s, zfsw ) 796 CALL wrk_dealloc( jpij, zf, dzf, zqns_ice_b, zerrit, zdifcase, zftrice, zihic, zghe ) 797 CALL wrk_dealloc( jpij,nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 ) 798 CALL wrk_dealloc( jpij,nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart = 0 ) 799 CALL wrk_dealloc( jpij,nlay_i+3, zindterm, zindtbis, zdiagbis ) 800 CALL wrk_dealloc( jpij,nlay_i+3,3, ztrid ) 801 CALL wrk_dealloc( jpij, zdq, zq_ini, zhfx_err ) 766 802 767 803 END SUBROUTINE lim_thd_dif … … 778 814 ! 779 815 INTEGER :: ji, jk ! dummy loop indices 780 REAL(wp) :: ztmelts , zindb! local scalar816 REAL(wp) :: ztmelts ! local scalar 781 817 !!------------------------------------------------------------------- 782 818 ! 783 819 DO jk = 1, nlay_i ! Sea ice energy of melting 784 820 DO ji = kideb, kiut 785 ztmelts = - tmut * s_i_b(ji,jk) + rtt 786 zindb = MAX( 0._wp , SIGN( 1._wp , -(t_i_b(ji,jk) - rtt) - epsi10 ) ) 787 q_i_b(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_b(ji,jk) ) & 788 & + lfus * ( 1.0 - zindb * ( ztmelts-rtt ) / MIN( t_i_b(ji,jk)-rtt, -epsi10 ) ) & 789 & - rcp * ( ztmelts-rtt ) ) 821 ztmelts = - tmut * s_i_1d(ji,jk) + rt0 822 t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts ) ! Force t_i_1d to be lower than melting point 823 ! (sometimes dif scheme produces abnormally high temperatures) 824 q_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) ) & 825 & + lfus * ( 1.0 - ( ztmelts-rt0 ) / ( t_i_1d(ji,jk) - rt0 ) ) & 826 & - rcp * ( ztmelts-rt0 ) ) 790 827 END DO 791 828 END DO 792 829 DO jk = 1, nlay_s ! Snow energy of melting 793 830 DO ji = kideb, kiut 794 q_s_ b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus )831 q_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) 795 832 END DO 796 833 END DO
Note: See TracChangeset
for help on using the changeset viewer.
