- Timestamp:
- 2017-09-13T18:46:56+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dif.F90
r8514 r8518 94 94 REAL(wp) :: ztsu_err = 1.e-5_wp ! range around which t_su is considered at 0C 95 95 REAL(wp) :: ztmelt_i ! ice melting temperature 96 REAL(wp) :: z hsu96 REAL(wp) :: z1_hsu 97 97 REAL(wp) :: zdti_max ! current maximal error on temperature 98 98 REAL(wp) :: zdti_bnd = 1.e-4_wp ! maximal authorized error on temperature 99 REAL(wp) :: zcpi ! Ice specific heat 100 REAL(wp) :: zhi ! 101 REAL(wp) :: zhfx_err, zdq ! diag errors on heat 99 102 100 103 REAL(wp), DIMENSION(jpij) :: isnow ! switch for presence (1) or absence (0) of snow 101 REAL(wp), DIMENSION(jpij) :: ztsub ! old surface temperature (before the iterative procedure ) 102 REAL(wp), DIMENSION(jpij) :: ztsubit ! surface temperature at previous iteration 104 REAL(wp), DIMENSION(jpij) :: ztsub ! surface temperature at previous iteration 103 105 REAL(wp), DIMENSION(jpij) :: zh_i ! ice layer thickness 104 106 REAL(wp), DIMENSION(jpij) :: zh_s ! snow layer thickness … … 106 108 REAL(wp), DIMENSION(jpij) :: zqns_ice_b ! solar radiation absorbed at the surface 107 109 REAL(wp), DIMENSION(jpij) :: zf ! surface flux function 108 REAL(wp), DIMENSION(jpij) :: dzf! derivative of the surface flux function110 REAL(wp), DIMENSION(jpij) :: zdqns_ice_b ! derivative of the surface flux function 109 111 REAL(wp), DIMENSION(jpij) :: zdti ! current error on temperature 110 REAL(wp), DIMENSION(jpij) :: zdifcase ! case of the equation resolution (1->4)111 112 REAL(wp), DIMENSION(jpij) :: zftrice ! solar radiation transmitted through the ice 112 REAL(wp), DIMENSION(jpij) :: zihic113 113 114 REAL(wp), DIMENSION(jpij,nlay_i) :: z_i ! Vertical cotes of the layers in the ice 115 REAL(wp), DIMENSION(jpij,nlay_s) :: z_s ! Vertical cotes of the layers in the snow 116 REAL(wp), DIMENSION(jpij,nlay_i) :: ztib ! Old temperature in the ice 117 REAL(wp), DIMENSION(jpij,nlay_s) :: ztsb ! Temporary temperature in the snow 118 REAL(wp), DIMENSION(jpij,nlay_i) :: ztitemp ! Temporary temperature in the ice to check the convergence 119 REAL(wp), DIMENSION(jpij,nlay_s) :: ztstemp ! Temporary temperature in the snow to check the convergence 114 120 REAL(wp), DIMENSION(jpij,0:nlay_i) :: ztcond_i ! Ice thermal conductivity 115 121 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradtr_i ! Radiation transmitted through the ice 116 122 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradab_i ! Radiation absorbed in the ice 117 123 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zkappa_i ! Kappa factor in the ice 118 REAL(wp), DIMENSION(jpij,0:nlay_i) :: ztib ! Old temperature in the ice119 124 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zeta_i ! Eta factor in the ice 120 REAL(wp), DIMENSION(jpij,0:nlay_i) :: ztitemp ! Temporary temperature in the ice to check the convergence121 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zspeche_i ! Ice specific heat122 REAL(wp), DIMENSION(jpij,0:nlay_i) :: z_i ! Vertical cotes of the layers in the ice123 125 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zradtr_s ! Radiation transmited through the snow 124 126 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zradab_s ! Radiation absorbed in the snow 125 127 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zkappa_s ! Kappa factor in the snow 126 128 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zeta_s ! Eta factor in the snow 127 REAL(wp), DIMENSION(jpij,0:nlay_s) :: ztstemp ! Temporary temperature in the snow to check the convergence128 REAL(wp), DIMENSION(jpij,0:nlay_s) :: ztsb ! Temporary temperature in the snow129 REAL(wp), DIMENSION(jpij,0:nlay_s) :: z_s ! Vertical cotes of the layers in the snow130 129 REAL(wp), DIMENSION(jpij,nlay_i+3) :: zindterm ! 'Ind'ependent term 131 130 REAL(wp), DIMENSION(jpij,nlay_i+3) :: zindtbis ! Temporary 'ind'ependent term 132 131 REAL(wp), DIMENSION(jpij,nlay_i+3) :: zdiagbis ! Temporary 'dia'gonal term 133 132 REAL(wp), DIMENSION(jpij,nlay_i+3,3) :: ztrid ! Tridiagonal system terms 134 REAL(wp), DIMENSION(jpij) :: z dq, zq_ini, zhfx_err! diag errors on heat133 REAL(wp), DIMENSION(jpij) :: zq_ini ! diag errors on heat 135 134 REAL(wp), DIMENSION(jpij) :: zghe ! G(he), th. conduct enhancement factor, mono-cat 136 135 … … 148 147 149 148 ! --- diag error on heat diffusion - PART 1 --- ! 150 zdq(:) = 0._wp ; zq_ini(:) = 0._wp151 149 DO ji = 1, nidx 152 150 zq_ini(ji) = ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i + & … … 167 165 ! Ice / snow layers 168 166 !-------------------- 169 170 z_s(:,0) = 0._wp ! vert. coord. of the up. lim. of the 1st snow layer 171 z_i(:,0) = 0._wp ! vert. coord. of the up. lim. of the 1st ice layer 172 173 DO jk = 1, nlay_s ! vert. coord of the up. lim. of the layer-th snow layer 174 DO ji = 1 , nidx 175 z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) * r1_nlay_s 176 END DO 177 END DO 178 179 DO jk = 1, nlay_i ! vert. coord of the up. lim. of the layer-th ice layer 180 DO ji = 1 , nidx 181 z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) * r1_nlay_i 167 z_s(1:nidx,1) = zh_s(1:nidx) 168 z_i(1:nidx,1) = zh_i(1:nidx) 169 DO jk = 2, nlay_s ! vert. coord of the up. lim. of the layer-th snow layer 170 DO ji = 1 , nidx 171 z_s(ji,jk) = z_s(ji,jk-1) + zh_s(ji) 172 END DO 173 END DO 174 175 DO jk = 2, nlay_i ! vert. coord of the up. lim. of the layer-th ice layer 176 DO ji = 1 , nidx 177 z_i(ji,jk) = z_i(ji,jk-1) + zh_i(ji) 182 178 END DO 183 179 END DO … … 187 183 !------------------------------------------------------------------------------| 188 184 ! 189 !------------------- 190 ! Computation of i0 191 !------------------- 192 ! i0 describes the fraction of solar radiation which does not contribute 193 ! to the surface energy budget but rather penetrates inside the ice. 194 ! We assume that no radiation is transmitted through the snow 195 ! If there is no no snow 196 ! zfsw = (1-i0).qsr_ice is absorbed at the surface 197 ! zftrice = io.qsr_ice is below the surface 198 ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice 199 ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover 200 zhsu = 0.1_wp ! threshold for the computation of i0 185 z1_hsu = 1._wp / 0.1_wp ! threshold for the computation of i0 201 186 DO ji = 1 , nidx 202 ! switches 203 isnow(ji) = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) ) 204 ! hs > 0, isnow = 1 205 zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) / zhsu ) ) 206 207 i0(ji) = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) ) 208 END DO 209 210 !------------------------------------------------------- 211 ! Solar radiation absorbed / transmitted at the surface 212 ! Derivative of the non solar flux 213 !------------------------------------------------------- 214 DO ji = 1 , nidx 215 zfsw (ji) = qsr_ice_1d(ji) * ( 1 - i0(ji) ) ! Shortwave radiation absorbed at surface 216 zftrice(ji) = qsr_ice_1d(ji) * i0(ji) ! Solar radiation transmitted below the surface layer 217 dzf (ji) = dqns_ice_1d(ji) ! derivative of incoming nonsolar flux 218 zqns_ice_b(ji) = qns_ice_1d(ji) ! store previous qns_ice_1d value 187 !------------------- 188 ! Computation of i0 189 !------------------- 190 ! i0 describes the fraction of solar radiation which does not contribute 191 ! to the surface energy budget but rather penetrates inside the ice. 192 ! We assume that no radiation is transmitted through the snow 193 ! If there is no no snow 194 ! zfsw = (1-i0).qsr_ice is absorbed at the surface 195 ! zftrice = io.qsr_ice is below the surface 196 ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice 197 ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover 198 zhi = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) * z1_hsu ) ) 199 i0(ji) = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zhi * fr2_i0_1d(ji) ) 200 201 !------------------------------------------------------- 202 ! Solar radiation absorbed / transmitted at the surface 203 ! Derivative of the non solar flux 204 !------------------------------------------------------- 205 zfsw (ji) = qsr_ice_1d(ji) * ( 1 - i0(ji) ) ! Shortwave radiation absorbed at surface 206 zftrice(ji) = qsr_ice_1d(ji) * i0(ji) ! Solar radiation transmitted below the surface layer 207 zdqns_ice_b(ji) = dqns_ice_1d(ji) ! derivative of incoming nonsolar flux 208 zqns_ice_b (ji) = qns_ice_1d(ji) ! store previous qns_ice_1d value 219 209 END DO 220 210 … … 222 212 ! Transmission - absorption of solar radiation in the ice 223 213 !--------------------------------------------------------- 224 225 DO ji = 1, nidx ! snow initialization 226 zradtr_s(ji,0) = zftrice(ji) ! radiation penetrating through snow 227 END DO 228 229 DO jk = 1, nlay_s ! Radiation through snow 214 zradtr_s(1:nidx,0) = zftrice(1:nidx) 215 DO jk = 1, nlay_s 230 216 DO ji = 1, nidx 231 217 ! ! radiation transmitted below the layer-th snow layer 232 zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,jk) )) )218 zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * z_s(ji,jk) ) 233 219 ! ! radiation absorbed by the layer-th snow layer 234 220 zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk) 235 221 END DO 236 222 END DO 237 238 DO ji = 1, nidx ! ice initialization 239 zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * isnow(ji) + zftrice(ji) * ( 1._wp - isnow(ji) ) 240 END DO 241 242 DO jk = 1, nlay_i ! Radiation through ice 223 224 zradtr_i(1:nidx,0) = zradtr_s(1:nidx,nlay_s) * isnow(1:nidx) + zftrice(1:nidx) * ( 1._wp - isnow(1:nidx) ) 225 DO jk = 1, nlay_i 243 226 DO ji = 1, nidx 244 227 ! ! radiation transmitted below the layer-th ice layer 245 zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * ( MAX ( 0._wp , z_i(ji,jk) )) )228 zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * z_i(ji,jk) ) 246 229 ! ! radiation absorbed by the layer-th ice layer 247 230 zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk) … … 249 232 END DO 250 233 251 DO ji = 1, nidx ! Radiation transmitted below the ice 252 ftr_ice_1d(ji) = zradtr_i(ji,nlay_i) 253 END DO 234 ftr_ice_1d(1:nidx) = zradtr_i(1:nidx,nlay_i) ! record radiation transmitted below the ice 254 235 255 236 !------------------------------------------------------------------------------| … … 257 238 !------------------------------------------------------------------------------| 258 239 ! 259 DO ji = 1, nidx ! Old surface temperature 260 ztsub (ji) = t_su_1d(ji) ! temperature at the beg of iter pr. 261 ztsubit(ji) = t_su_1d(ji) ! temperature at the previous iter 262 t_su_1d(ji) = MIN( t_su_1d(ji), rt0 - ztsu_err ) ! necessary 263 zdti (ji) = 1000._wp ! initial value of error 264 END DO 265 266 DO jk = 1, nlay_s ! Old snow temperature 267 DO ji = 1 , nidx 268 ztsb(ji,jk) = t_s_1d(ji,jk) 269 END DO 270 END DO 271 272 DO jk = 1, nlay_i ! Old ice temperature 273 DO ji = 1 , nidx 274 ztib(ji,jk) = t_i_1d(ji,jk) 275 END DO 276 END DO 240 ztsub (1:nidx) = t_su_1d(1:nidx) ! temperature at the previous iter 241 t_su_1d(1:nidx) = MIN( t_su_1d(1:nidx), rt0 - ztsu_err ) ! necessary 242 zdti (1:nidx) = 1000._wp ! initial value of error 243 ztsb (1:nidx,:) = t_s_1d(1:nidx,:) ! Old snow temperature 244 ztib (1:nidx,:) = t_i_1d(1:nidx,:) ! Old ice temperature 277 245 278 246 iconv = 0 ! number of iterations … … 289 257 IF( ln_cndi_U64 ) THEN !-- Untersteiner (1964) formula: k = k0 + beta.S/T 290 258 DO ji = 1 , nidx 291 ztcond_i(ji,0) = rcdic + zbeta * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) 292 ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin ) 259 ztcond_i(ji,0) = MAX( zkimin, rcdic + zbeta * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) ) 293 260 END DO 294 261 DO jk = 1, nlay_i-1 295 262 DO ji = 1 , nidx 296 ztcond_i(ji,jk) = rcdic + zbeta * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) / & 297 MIN(-2.0_wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0) 298 ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin ) 263 ztcond_i(ji,jk) = MAX( zkimin, rcdic + zbeta * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) / & 264 & MIN(-2.0_wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0) ) 299 265 END DO 300 266 END DO … … 302 268 ELSEIF( ln_cndi_P07 ) THEN !-- Pringle et al formula: k = k0 + beta1.S/T - beta2.T 303 269 DO ji = 1 , nidx 304 ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) & 305 & - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) 306 ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin ) 270 ztcond_i(ji,0) = MAX( zkimin, rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) & 271 & - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) ) 307 272 END DO 308 273 DO jk = 1, nlay_i-1 309 274 DO ji = 1 , nidx 310 ztcond_i(ji,jk) = rcdic +&275 ztcond_i(ji,jk) = MAX( zkimin, rcdic + & 311 276 & 0.09_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) & 312 277 & / MIN( -2._wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0 ) & 313 & - 0.0055_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0 * rt0 ) 314 ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin ) 278 & - 0.0055_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0 * rt0 ) ) 315 279 END DO 316 280 END DO 317 281 DO ji = 1 , nidx 318 ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) & 319 & - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 320 ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin ) 282 ztcond_i(ji,nlay_i) = MAX( zkimin, rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) & 283 & - 0.011_wp * ( t_bo_1d(ji) - rt0 ) ) 321 284 END DO 322 285 ENDIF … … 335 298 SELECT CASE ( nn_monocat ) 336 299 337 CASE (1,3) ! LIM3300 CASE (1,3) 338 301 339 302 zepsilon = 0.1_wp … … 404 367 DO jk = 1, nlay_i 405 368 DO ji = 1 , nidx 406 ztitemp(ji,jk) 407 z speche_i(ji,jk)= cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztib(ji,jk) - rt0 ), epsi10 )408 zeta_i(ji,jk) = rdt_ice / MAX( rhoic * zspeche_i(ji,jk)* zh_i(ji), epsi10 )369 ztitemp(ji,jk) = t_i_1d(ji,jk) 370 zcpi = cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztib(ji,jk) - rt0 ), epsi10 ) 371 zeta_i(ji,jk) = rdt_ice / MAX( rhoic * zcpi * zh_i(ji), epsi10 ) 409 372 END DO 410 373 END DO … … 425 388 DO ji = 1 , nidx 426 389 ! update of the non solar flux according to the update in T_su 427 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub it(ji) )390 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) 428 391 END DO 429 392 ENDIF … … 507 470 ! case 1 : no surface melting - snow present | 508 471 !------------------------------------------------------------------------------| 509 zdifcase(ji) = 1.0510 472 numeqmin(ji) = 1 511 473 numeqmax(ji) = nlay_i + nlay_s + 1 … … 513 475 !!surface equation 514 476 ztrid(ji,1,1) = 0.0 515 ztrid(ji,1,2) = dzf(ji) - zg1s * zkappa_s(ji,0)477 ztrid(ji,1,2) = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0) 516 478 ztrid(ji,1,3) = zg1s * zkappa_s(ji,0) 517 zindterm(ji,1) = dzf(ji) * t_su_1d(ji) - zf(ji)479 zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zf(ji) 518 480 519 481 !!first layer of snow equation … … 529 491 !------------------------------------------------------------------------------| 530 492 ! 531 zdifcase(ji) = 2.0532 493 numeqmin(ji) = 2 533 494 numeqmax(ji) = nlay_i + nlay_s + 1 … … 552 513 !------------------------------------------------------------------------------| 553 514 ! 554 zdifcase(ji) = 3.0555 515 numeqmin(ji) = nlay_s + 1 556 516 numeqmax(ji) = nlay_i + nlay_s + 1 … … 558 518 !!surface equation 559 519 ztrid(ji,numeqmin(ji),1) = 0.0 560 ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0)*zg1520 ztrid(ji,numeqmin(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0)*zg1 561 521 ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0)*zg1 562 zindterm(ji,numeqmin(ji)) = dzf(ji)*t_su_1d(ji) - zf(ji)522 zindterm(ji,numeqmin(ji)) = zdqns_ice_b(ji)*t_su_1d(ji) - zf(ji) 563 523 564 524 !!first layer of ice equation … … 572 532 IF ( nlay_i == 1 ) THEN 573 533 ztrid(ji,numeqmin(ji),1) = 0.0 574 ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0) * 2.0534 ztrid(ji,numeqmin(ji),2) = zdqns_ice_b(ji) - zkappa_i(ji,0) * 2.0 575 535 ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0) * 2.0 576 536 ztrid(ji,numeqmin(ji)+1,1) = -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1) … … 589 549 !------------------------------------------------------------------------------| 590 550 ! 591 zdifcase(ji) = 4.0592 551 numeqmin(ji) = nlay_s + 2 593 552 numeqmax(ji) = nlay_i + nlay_s + 1 … … 662 621 ! surface temperature 663 622 isnow(ji) = 1._wp - MAX( 0._wp , SIGN( 1._wp , -ht_s_1d(ji) ) ) 664 ztsub it(ji) = t_su_1d(ji)623 ztsub(ji) = t_su_1d(ji) 665 624 IF( t_su_1d(ji) < rt0 ) & 666 625 t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3) * & … … 676 635 DO ji = 1 , nidx 677 636 t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , 190._wp ) 678 zdti (ji) = ABS( t_su_1d(ji) - ztsub it(ji) )637 zdti (ji) = ABS( t_su_1d(ji) - ztsub(ji) ) 679 638 END DO 680 639 … … 754 713 755 714 ! --- diag conservation imbalance on heat diffusion - PART 2 --- ! 715 ! hfx_dif = Heat flux used to warm/cool ice in W.m-2 716 ! zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation 756 717 DO ji = 1, nidx 757 zdq(ji) = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i + & 758 & SUM( e_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s ) 718 zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i + & 719 & SUM( e_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s ) 720 759 721 IF( t_su_1d(ji) < rt0 ) THEN ! case T_su < 0degC 760 zhfx_err (ji) = qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice722 zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) 761 723 ELSE ! case T_su = 0degC 762 zhfx_err (ji) = fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice724 zhfx_err = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) 763 725 ENDIF 764 hfx_ err_1d(ji) = hfx_err_1d(ji) + zhfx_err(ji)* a_i_1d(ji)726 hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji) 765 727 766 728 ! total heat that is sent to the ocean (i.e. not used in the heat diffusion equation) 767 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err(ji) * a_i_1d(ji) 768 END DO 769 770 !----------------------------------------- 771 ! Heat flux used to warm/cool ice in W.m-2 772 !----------------------------------------- 773 DO ji = 1, nidx 774 IF( t_su_1d(ji) < rt0 ) THEN ! case T_su < 0degC 775 hfx_dif_1d(ji) = hfx_dif_1d(ji) + & 776 & ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji) 777 ELSE ! case T_su = 0degC 778 hfx_dif_1d(ji) = hfx_dif_1d(ji) + & 779 & ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji) 780 ENDIF 781 ! correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation 782 hfx_dif_1d(ji) = hfx_dif_1d(ji) - zhfx_err(ji) * a_i_1d(ji) 729 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err 730 783 731 END DO 784 732 !
Note: See TracChangeset
for help on using the changeset viewer.