Changeset 921 for trunk/NEMO/LIM_SRC_3/limthd_dif.F90
- Timestamp:
- 2008-05-13T10:28:52+02:00 (16 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/LIM_SRC_3/limthd_dif.F90 (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/LIM_SRC_3/limthd_dif.F90
r869 r921 22 22 USE par_ice 23 23 USE lib_mpp 24 24 25 25 IMPLICIT NONE 26 26 PRIVATE … … 44 44 45 45 SUBROUTINE lim_thd_dif( kideb , kiut , jl ) 46 !!------------------------------------------------------------------ 47 !! *** ROUTINE lim_thd_dif *** 48 !! ** Purpose : 49 !! This routine determines the time evolution of snow and sea-ice 50 !! temperature profiles. 51 !! ** Method : 52 !! This is done by solving the heat equation diffusion with 53 !! a Neumann boundary condition at the surface and a Dirichlet one 54 !! at the bottom. Solar radiation is partially absorbed into the ice. 55 !! The specific heat and thermal conductivities depend on ice salinity 56 !! and temperature to take into account brine pocket melting. The 57 !! numerical 58 !! scheme is an iterative Crank-Nicolson on a non-uniform multilayer grid 59 !! in the ice and snow system. 60 !! 61 !! The successive steps of this routine are 62 !! 1. Thermal conductivity at the interfaces of the ice layers 63 !! 2. Internal absorbed radiation 64 !! 3. Scale factors due to non-uniform grid 65 !! 4. Kappa factors 66 !! Then iterative procedure begins 67 !! 5. specific heat in the ice 68 !! 6. eta factors 69 !! 7. surface flux computation 70 !! 8. tridiagonal system terms 71 !! 9. solving the tridiagonal system with Gauss elimination 72 !! Iterative procedure ends according to a criterion on evolution 73 !! of temperature 74 !! 75 !! ** Arguments : 76 !! kideb , kiut : Starting and ending points on which the 77 !! the computation is applied 78 !! 79 !! ** Inputs / Ouputs : (global commons) 80 !! surface temperature : t_su_b 81 !! ice/snow temperatures : t_i_b, t_s_b 82 !! ice salinities : s_i_b 83 !! number of layers in the ice/snow: nlay_i, nlay_s 84 !! profile of the ice/snow layers : z_i, z_s 85 !! total ice/snow thickness : ht_i_b, ht_s_b 86 !! 87 !! ** External : 88 !! 89 !! ** References : 90 !! 91 !! ** History : 92 !! (02-2003) Martin Vancoppenolle, Louvain-la-Neuve, Belgium 93 !! (06-2005) Martin Vancoppenolle, 3d version 94 !! (11-2006) Vectorized by Xavier Fettweis (UCL-ASTR) 95 !! (04-2007) Energy conservation tested by M. Vancoppenolle 96 !! 97 !!------------------------------------------------------------------ 98 !! * Arguments 99 100 INTEGER , INTENT (in) :: & 101 kideb , & ! Start point on which the the computation is applied 102 kiut , & ! End point on which the the computation is applied 103 jl ! Category number 104 105 !! * Local variables 106 INTEGER :: ji, & ! spatial loop index 107 zji, zjj, & ! temporary dummy loop index 108 numeq, & ! current reference number of equation 109 layer, & ! vertical dummy loop index 110 nconv, & ! number of iterations in iterative procedure 111 minnumeqmin, & ! 112 maxnumeqmax 113 114 INTEGER , DIMENSION(jpij) :: & 115 numeqmin, & ! reference number of top equation 116 numeqmax, & ! reference number of bottom equation 117 isnow ! switch for presence (1) or absence (0) of snow 118 119 !! * New local variables 120 REAL(wp) , DIMENSION(jpij,0:nlay_i) :: & 121 ztcond_i, & !Ice thermal conductivity 122 zradtr_i, & !Radiation transmitted through the ice 123 zradab_i, & !Radiation absorbed in the ice 124 zkappa_i !Kappa factor in the ice 125 126 REAL(wp) , DIMENSION(jpij,0:nlay_s) :: & 127 zradtr_s, & !Radiation transmited through the snow 128 zradab_s, & !Radiation absorbed in the snow 129 zkappa_s !Kappa factor in the snow 130 131 REAL(wp) , DIMENSION(jpij,0:nlay_i) :: & 132 ztiold, & !Old temperature in the ice 133 zeta_i, & !Eta factor in the ice 134 ztitemp, & !Temporary temperature in the ice to check the convergence 135 zspeche_i, & !Ice specific heat 136 z_i !Vertical cotes of the layers in the ice 137 138 REAL(wp) , DIMENSION(jpij,0:nlay_s) :: & 139 zeta_s, & !Eta factor in the snow 140 ztstemp, & !Temporary temperature in the snow to check the convergence 141 ztsold, & !Temporary temperature in the snow 142 z_s !Vertical cotes of the layers in the snow 143 144 REAL(wp) , DIMENSION(jpij,jkmax+2) :: & 145 zindterm, & ! Independent term 146 zindtbis, & ! temporary independent term 147 zdiagbis 148 149 REAL(wp) , DIMENSION(jpij,jkmax+2,3) :: & 150 ztrid ! tridiagonal system terms 151 152 REAL(wp), DIMENSION(jpij) :: & 153 ztfs , & ! ice melting point 154 ztsuold , & ! old surface temperature (before the iterative 155 ! procedure ) 156 ztsuoldit, & ! surface temperature at previous iteration 157 zh_i , & !ice layer thickness 158 zh_s , & !snow layer thickness 159 zfsw , & !solar radiation absorbed at the surface 160 zf , & ! surface flux function 161 dzf ! derivative of the surface flux function 162 163 REAL(wp) :: & ! constant values 164 zeps = 1.0e-10, & ! 165 zg1s = 2.0, & !: for the tridiagonal system 166 zg1 = 2.0, & 167 zgamma = 18009.0, & !: for specific heat 168 zbeta = 0.117, & !: for thermal conductivity (could be 0.13) 169 zraext_s = 1.0e08, & !: extinction coefficient of radiation in the snow 170 zkimin = 0.10 , & !: minimum ice thermal conductivity 171 zht_smin = 1.0e-4 !: minimum snow depth 172 173 REAL(wp) :: & ! local variables 174 ztmelt_i, & ! ice melting temperature 175 zerritmax ! current maximal error on temperature 176 177 REAL(wp), DIMENSION(jpij) :: & 178 zerrit, & ! current error on temperature 179 zdifcase, & ! case of the equation resolution (1->4) 180 zftrice, & ! solar radiation transmitted through the ice 181 zihic, zhsu 182 183 !!-- End of declarations 184 !!---------------------------------------------------------------------------------------------- 185 186 IF(lwp) WRITE(numout,*)'lim_thd_dif : Heat diffusion in sea ice for cat :', jl 187 188 ! 189 !------------------------------------------------------------------------------! 190 ! 1) Initialization ! 191 !------------------------------------------------------------------------------! 192 ! 193 DO ji = kideb , kiut 194 ! is there snow or not 195 isnow(ji)= INT ( 1.0 - MAX( 0.0 , SIGN (1.0, - ht_s_b(ji) ) ) ) 196 ! surface temperature of fusion 197 ztfs(ji) = isnow(ji) * rtt + (1.0-isnow(ji)) * rtt 198 ! layer thickness 199 zh_i(ji) = ht_i_b(ji) / nlay_i 200 zh_s(ji) = ht_s_b(ji) / nlay_s 201 END DO 202 203 !-------------------- 204 ! Ice / snow layers 205 !-------------------- 206 207 z_s(:,0) = 0.0 ! vert. coord. of the up. lim. of the 1st snow layer 208 z_i(:,0) = 0.0 ! vert. coord. of the up. lim. of the 1st ice layer 209 210 DO layer = 1, nlay_s 211 DO ji = kideb , kiut 212 ! vert. coord of the up. lim. of the layer-th snow layer 213 z_s(ji,layer) = z_s(ji,layer-1) + ht_s_b(ji) / nlay_s 214 END DO 215 END DO 216 217 DO layer = 1, nlay_i 218 DO ji = kideb , kiut 219 ! vert. coord of the up. lim. of the layer-th ice layer 220 z_i(ji,layer) = z_i(ji,layer-1) + ht_i_b(ji) / nlay_i 221 END DO 222 END DO 223 ! 224 !------------------------------------------------------------------------------| 225 ! 2) Radiations | 226 !------------------------------------------------------------------------------| 227 ! 228 !------------------- 229 ! Computation of i0 230 !------------------- 231 ! i0 describes the fraction of solar radiation which does not contribute 232 ! to the surface energy budget but rather penetrates inside the ice. 233 ! We assume that no radiation is transmitted through the snow 234 ! If there is no no snow 235 ! zfsw = (1-i0).qsr_ice is absorbed at the surface 236 ! zftrice = io.qsr_ice is below the surface 237 ! fstbif = io.qsr_ice.exp(-k(h_i)) transmitted below the ice 238 239 DO ji = kideb , kiut 240 ! switches 241 isnow(ji) = INT ( 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_s_b(ji) ) ) ) 242 ! hs > 0, isnow = 1 243 zhsu(ji) = hnzst !threshold for the computation of i0 244 zihic(ji) = MAX( zzero , 1.0 - ( ht_i_b(ji) / zhsu(ji) ) ) 245 246 i0(ji) = ( 1.0 - isnow(ji) ) * & 247 ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) ) 248 !fr1_i0_1d = i0 for a thin ice surface 249 !fr1_i0_2d = i0 for a thick ice surface 250 ! a function of the cloud cover 251 ! 252 !i0(ji) = (1.0-FLOAT(isnow(ji)))*3.0/(100*ht_s_b(ji)+10.0) 253 !formula used in Cice 254 END DO 255 256 !------------------------------------------------------- 257 ! Solar radiation absorbed / transmitted at the surface 258 ! Derivative of the non solar flux 259 !------------------------------------------------------- 260 DO ji = kideb , kiut 261 262 ! Shortwave radiation absorbed at surface 263 zfsw(ji) = qsr_ice_1d(ji) * ( 1 - i0(ji) ) 264 265 ! Solar radiation transmitted below the surface layer 266 zftrice(ji)= qsr_ice_1d(ji) * i0(ji) 267 268 ! derivative of incoming nonsolar flux 269 dzf(ji) = dqns_ice_1d(ji) 270 271 END DO 272 273 !--------------------------------------------------------- 274 ! Transmission - absorption of solar radiation in the ice 275 !--------------------------------------------------------- 276 277 DO ji = kideb , kiut 278 ! Initialization 279 zradtr_s(ji,0) = zftrice(ji) ! radiation penetrating through snow 280 END DO 281 282 ! Radiation through snow 283 DO layer = 1, nlay_s 284 DO ji = kideb , kiut 285 ! radiation transmitted below the layer-th snow layer 286 zradtr_s(ji,layer) = zradtr_s(ji,0) * EXP ( - zraext_s * ( MAX ( 0.0 , & 287 z_s(ji,layer) ) ) ) 288 ! radiation absorbed by the layer-th snow layer 289 zradab_s(ji,layer) = zradtr_s(ji,layer-1) - zradtr_s(ji,layer) 290 END DO 291 END DO 292 293 ! Radiation through ice 294 DO ji = kideb , kiut 295 zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * isnow(ji) + & 296 zftrice(ji) * ( 1 - isnow(ji) ) 297 END DO 298 299 DO layer = 1, nlay_i 300 DO ji = kideb , kiut 301 ! radiation transmitted below the layer-th ice layer 302 zradtr_i(ji,layer) = zradtr_i(ji,0) * EXP ( - kappa_i * ( MAX ( 0.0 , & 303 z_i(ji,layer) ) ) ) 304 ! radiation absorbed by the layer-th ice layer 305 zradab_i(ji,layer) = zradtr_i(ji,layer-1) - zradtr_i(ji,layer) 306 END DO 307 END DO 308 309 ! Radiation transmitted below the ice 310 DO ji = kideb , kiut 311 fstbif_1d(ji) = fstbif_1d(ji) + & 312 zradtr_i(ji,nlay_i) * a_i_b(ji) / at_i_b(ji) 313 END DO 314 315 ! +++++ 316 ! just to check energy conservation 317 DO ji = kideb , kiut 318 zji = MOD( npb(ji) - 1, jpi ) + 1 319 zjj = ( npb(ji) - 1 ) / jpi + 1 320 fstroc(zji,zjj,jl) = & 321 zradtr_i(ji,nlay_i) 322 END DO 323 ! +++++ 324 325 DO layer = 1, nlay_i 326 DO ji = kideb , kiut 327 radab(ji,layer) = zradab_i(ji,layer) 328 END DO 329 END DO 330 331 332 ! 333 !------------------------------------------------------------------------------| 334 ! 3) Iterative procedure begins | 335 !------------------------------------------------------------------------------| 336 ! 337 ! Old surface temperature 338 DO ji = kideb, kiut 339 ztsuold(ji) = t_su_b(ji) ! temperature at the beg of iter pr. 340 ztsuoldit(ji) = t_su_b(ji) ! temperature at the previous iter 341 t_su_b(ji) = MIN(t_su_b(ji),ztfs(ji)-0.00001) !necessary 342 zerrit(ji) = 1000.0 ! initial value of error 343 END DO 344 !RB Min global ?? 345 346 ! Old snow temperature 347 DO layer = 1, nlay_s 348 DO ji = kideb , kiut 349 ztsold(ji,layer) = t_s_b(ji,layer) 350 END DO 351 END DO 352 353 ! Old ice temperature 354 DO layer = 1, nlay_i 355 DO ji = kideb , kiut 356 ztiold(ji,layer) = t_i_b(ji,layer) 357 END DO 358 END DO 359 360 nconv = 0 ! number of iterations 361 zerritmax = 1000.0 ! maximal value of error on all points 362 363 DO WHILE ((zerritmax > maxer_i_thd).AND.(nconv < nconv_i_thd)) 364 365 nconv = nconv+1 366 367 ! 368 !------------------------------------------------------------------------------| 369 ! 4) Sea ice thermal conductivity | 370 !------------------------------------------------------------------------------| 371 ! 372 IF ( thcon_i_swi .EQ. 0 ) THEN 373 ! Untersteiner (1964) formula 374 DO ji = kideb , kiut 375 ztcond_i(ji,0) = rcdic + zbeta*s_i_b(ji,1) / & 376 MIN(-zeps,t_i_b(ji,1)-rtt) 377 ztcond_i(ji,0) = MAX(ztcond_i(ji,0),zkimin) 378 END DO 379 ENDIF 380 381 IF ( thcon_i_swi .EQ. 1 ) THEN 382 ! Pringle et al formula included, 383 ! 2.11 + 0.09 S/T - 0.011.T 384 DO ji = kideb , kiut 385 ztcond_i(ji,0) = rcdic + 0.09*s_i_b(ji,1) / & 386 MIN(-zeps,t_i_b(ji,1)-rtt) - & 387 0.011* ( t_i_b(ji,1) - rtt ) 388 ztcond_i(ji,0) = MAX(ztcond_i(ji,0),zkimin) 389 END DO 390 ENDIF 391 392 IF ( thcon_i_swi .EQ. 0 ) THEN ! Untersteiner 393 DO layer = 1, nlay_i-1 394 DO ji = kideb , kiut 395 ztcond_i(ji,layer) = rcdic + zbeta*( s_i_b(ji,layer) & 396 + s_i_b(ji,layer+1) ) / MIN(-zeps, & 397 t_i_b(ji,layer)+t_i_b(ji,layer+1)-2.0*rtt) 398 ztcond_i(ji,layer) = MAX(ztcond_i(ji,layer),zkimin) 399 END DO 400 END DO 401 ENDIF 402 403 IF ( thcon_i_swi .EQ. 1 ) THEN ! Pringle 404 DO layer = 1, nlay_i-1 405 DO ji = kideb , kiut 406 ztcond_i(ji,layer) = rcdic + 0.09*( s_i_b(ji,layer) & 407 + s_i_b(ji,layer+1) ) / MIN(-zeps, & 408 t_i_b(ji,layer)+t_i_b(ji,layer+1)-2.0*rtt) - & 409 0.011* ( t_i_b(ji,layer) + t_i_b(ji,layer+1) - 2.0*rtt ) 410 ztcond_i(ji,layer) = MAX(ztcond_i(ji,layer),zkimin) 411 END DO 412 END DO 413 ENDIF 414 415 IF ( thcon_i_swi .EQ. 0 ) THEN ! Untersteiner 416 DO ji = kideb , kiut 417 ztcond_i(ji,nlay_i) = rcdic + zbeta*s_i_b(ji,nlay_i) / & 418 MIN(-zeps,t_bo_b(ji)-rtt) 419 ztcond_i(ji,nlay_i) = MAX(ztcond_i(ji,nlay_i),zkimin) 420 END DO 421 ENDIF 422 423 IF ( thcon_i_swi .EQ. 1 ) THEN ! Pringle 424 DO ji = kideb , kiut 425 ztcond_i(ji,nlay_i) = rcdic + 0.09*s_i_b(ji,nlay_i) / & 426 MIN(-zeps,t_bo_b(ji)-rtt) - & 427 0.011* ( t_bo_b(ji) - rtt ) 428 ztcond_i(ji,nlay_i) = MAX(ztcond_i(ji,nlay_i),zkimin) 429 END DO 430 ENDIF 431 ! 432 !------------------------------------------------------------------------------| 433 ! 5) kappa factors | 434 !------------------------------------------------------------------------------| 435 ! 436 DO ji = kideb, kiut 437 438 !-- Snow kappa factors 439 zkappa_s(ji,0) = rcdsn / MAX(zeps,zh_s(ji)) 440 zkappa_s(ji,nlay_s) = rcdsn / MAX(zeps,zh_s(ji)) 441 END DO 442 443 DO layer = 1, nlay_s-1 444 DO ji = kideb , kiut 445 zkappa_s(ji,layer) = 2.0 * rcdsn / & 446 MAX(zeps,2.0*zh_s(ji)) 447 END DO 448 END DO 449 450 DO layer = 1, nlay_i-1 451 DO ji = kideb , kiut 452 !-- Ice kappa factors 453 zkappa_i(ji,layer) = 2.0*ztcond_i(ji,layer)/ & 454 MAX(zeps,2.0*zh_i(ji)) 455 END DO 456 END DO 457 458 DO ji = kideb , kiut 459 zkappa_i(ji,0) = ztcond_i(ji,0)/MAX(zeps,zh_i(ji)) 460 zkappa_i(ji,nlay_i) = ztcond_i(ji,nlay_i) / MAX(zeps,zh_i(ji)) 461 !-- Interface 462 zkappa_s(ji,nlay_s) = 2.0*rcdsn*ztcond_i(ji,0)/MAX(zeps, & 463 (ztcond_i(ji,0)*zh_s(ji) + rcdsn*zh_i(ji))) 464 zkappa_i(ji,0) = zkappa_s(ji,nlay_s)*isnow(ji) & 465 + zkappa_i(ji,0)*(1.0-isnow(ji)) 466 END DO 467 ! 468 !------------------------------------------------------------------------------| 469 ! 6) Sea ice specific heat, eta factors | 470 !------------------------------------------------------------------------------| 471 ! 472 DO layer = 1, nlay_i 473 DO ji = kideb , kiut 474 ztitemp(ji,layer) = t_i_b(ji,layer) 475 zspeche_i(ji,layer) = cpic + zgamma*s_i_b(ji,layer)/ & 476 MAX((t_i_b(ji,layer)-rtt)*(ztiold(ji,layer)-rtt),zeps) 477 zeta_i(ji,layer) = rdt_ice / MAX(rhoic*zspeche_i(ji,layer)*zh_i(ji), & 478 zeps) 479 END DO 480 END DO 481 482 DO layer = 1, nlay_s 483 DO ji = kideb , kiut 484 ztstemp(ji,layer) = t_s_b(ji,layer) 485 zeta_s(ji,layer) = rdt_ice / MAX(rhosn*cpic*zh_s(ji),zeps) 486 END DO 487 END DO 488 ! 489 !------------------------------------------------------------------------------| 490 ! 7) surface flux computation | 491 !------------------------------------------------------------------------------| 492 ! 493 DO ji = kideb , kiut 494 495 ! update of the non solar flux according to the update in T_su 496 qnsr_ice_1d(ji) = qnsr_ice_1d(ji) + dqns_ice_1d(ji) * & 497 ( t_su_b(ji) - ztsuoldit(ji) ) 498 499 ! update incoming flux 500 zf(ji) = zfsw(ji) & ! net absorbed solar radiation 501 + qnsr_ice_1d(ji) ! non solar total flux 502 ! (LWup, LWdw, SH, LH) 503 504 END DO 505 506 ! 507 !------------------------------------------------------------------------------| 508 ! 8) tridiagonal system terms | 509 !------------------------------------------------------------------------------| 510 ! 511 !!layer denotes the number of the layer in the snow or in the ice 512 !!numeq denotes the reference number of the equation in the tridiagonal 513 !!system, terms of tridiagonal system are indexed as following : 514 !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one 515 516 !!ice interior terms (top equation has the same form as the others) 517 518 DO numeq=1,jkmax+2 519 DO ji = kideb , kiut 520 ztrid(ji,numeq,1) = 0. 521 ztrid(ji,numeq,2) = 0. 522 ztrid(ji,numeq,3) = 0. 523 zindterm(ji,numeq)= 0. 524 zindtbis(ji,numeq)= 0. 525 zdiagbis(ji,numeq)= 0. 526 ENDDO 527 ENDDO 528 529 DO numeq = nlay_s + 2, nlay_s + nlay_i 530 DO ji = kideb , kiut 531 layer = numeq - nlay_s - 1 532 ztrid(ji,numeq,1) = - zeta_i(ji,layer)*zkappa_i(ji,layer-1) 533 ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,layer)*(zkappa_i(ji,layer-1) + & 534 zkappa_i(ji,layer)) 535 ztrid(ji,numeq,3) = - zeta_i(ji,layer)*zkappa_i(ji,layer) 536 zindterm(ji,numeq) = ztiold(ji,layer) + zeta_i(ji,layer)* & 537 zradab_i(ji,layer) 538 END DO 539 ENDDO 540 541 numeq = nlay_s + nlay_i + 1 542 DO ji = kideb , kiut 46 !!------------------------------------------------------------------ 47 !! *** ROUTINE lim_thd_dif *** 48 !! ** Purpose : 49 !! This routine determines the time evolution of snow and sea-ice 50 !! temperature profiles. 51 !! ** Method : 52 !! This is done by solving the heat equation diffusion with 53 !! a Neumann boundary condition at the surface and a Dirichlet one 54 !! at the bottom. Solar radiation is partially absorbed into the ice. 55 !! The specific heat and thermal conductivities depend on ice salinity 56 !! and temperature to take into account brine pocket melting. The 57 !! numerical 58 !! scheme is an iterative Crank-Nicolson on a non-uniform multilayer grid 59 !! in the ice and snow system. 60 !! 61 !! The successive steps of this routine are 62 !! 1. Thermal conductivity at the interfaces of the ice layers 63 !! 2. Internal absorbed radiation 64 !! 3. Scale factors due to non-uniform grid 65 !! 4. Kappa factors 66 !! Then iterative procedure begins 67 !! 5. specific heat in the ice 68 !! 6. eta factors 69 !! 7. surface flux computation 70 !! 8. tridiagonal system terms 71 !! 9. solving the tridiagonal system with Gauss elimination 72 !! Iterative procedure ends according to a criterion on evolution 73 !! of temperature 74 !! 75 !! ** Arguments : 76 !! kideb , kiut : Starting and ending points on which the 77 !! the computation is applied 78 !! 79 !! ** Inputs / Ouputs : (global commons) 80 !! surface temperature : t_su_b 81 !! ice/snow temperatures : t_i_b, t_s_b 82 !! ice salinities : s_i_b 83 !! number of layers in the ice/snow: nlay_i, nlay_s 84 !! profile of the ice/snow layers : z_i, z_s 85 !! total ice/snow thickness : ht_i_b, ht_s_b 86 !! 87 !! ** External : 88 !! 89 !! ** References : 90 !! 91 !! ** History : 92 !! (02-2003) Martin Vancoppenolle, Louvain-la-Neuve, Belgium 93 !! (06-2005) Martin Vancoppenolle, 3d version 94 !! (11-2006) Vectorized by Xavier Fettweis (UCL-ASTR) 95 !! (04-2007) Energy conservation tested by M. Vancoppenolle 96 !! 97 !!------------------------------------------------------------------ 98 !! * Arguments 99 100 INTEGER , INTENT (in) :: & 101 kideb , & ! Start point on which the the computation is applied 102 kiut , & ! End point on which the the computation is applied 103 jl ! Category number 104 105 !! * Local variables 106 INTEGER :: ji, & ! spatial loop index 107 zji, zjj, & ! temporary dummy loop index 108 numeq, & ! current reference number of equation 109 layer, & ! vertical dummy loop index 110 nconv, & ! number of iterations in iterative procedure 111 minnumeqmin, & ! 112 maxnumeqmax 113 114 INTEGER , DIMENSION(jpij) :: & 115 numeqmin, & ! reference number of top equation 116 numeqmax, & ! reference number of bottom equation 117 isnow ! switch for presence (1) or absence (0) of snow 118 119 !! * New local variables 120 REAL(wp) , DIMENSION(jpij,0:nlay_i) :: & 121 ztcond_i, & !Ice thermal conductivity 122 zradtr_i, & !Radiation transmitted through the ice 123 zradab_i, & !Radiation absorbed in the ice 124 zkappa_i !Kappa factor in the ice 125 126 REAL(wp) , DIMENSION(jpij,0:nlay_s) :: & 127 zradtr_s, & !Radiation transmited through the snow 128 zradab_s, & !Radiation absorbed in the snow 129 zkappa_s !Kappa factor in the snow 130 131 REAL(wp) , DIMENSION(jpij,0:nlay_i) :: & 132 ztiold, & !Old temperature in the ice 133 zeta_i, & !Eta factor in the ice 134 ztitemp, & !Temporary temperature in the ice to check the convergence 135 zspeche_i, & !Ice specific heat 136 z_i !Vertical cotes of the layers in the ice 137 138 REAL(wp) , DIMENSION(jpij,0:nlay_s) :: & 139 zeta_s, & !Eta factor in the snow 140 ztstemp, & !Temporary temperature in the snow to check the convergence 141 ztsold, & !Temporary temperature in the snow 142 z_s !Vertical cotes of the layers in the snow 143 144 REAL(wp) , DIMENSION(jpij,jkmax+2) :: & 145 zindterm, & ! Independent term 146 zindtbis, & ! temporary independent term 147 zdiagbis 148 149 REAL(wp) , DIMENSION(jpij,jkmax+2,3) :: & 150 ztrid ! tridiagonal system terms 151 152 REAL(wp), DIMENSION(jpij) :: & 153 ztfs , & ! ice melting point 154 ztsuold , & ! old surface temperature (before the iterative 155 ! procedure ) 156 ztsuoldit, & ! surface temperature at previous iteration 157 zh_i , & !ice layer thickness 158 zh_s , & !snow layer thickness 159 zfsw , & !solar radiation absorbed at the surface 160 zf , & ! surface flux function 161 dzf ! derivative of the surface flux function 162 163 REAL(wp) :: & ! constant values 164 zeps = 1.0e-10, & ! 165 zg1s = 2.0, & !: for the tridiagonal system 166 zg1 = 2.0, & 167 zgamma = 18009.0, & !: for specific heat 168 zbeta = 0.117, & !: for thermal conductivity (could be 0.13) 169 zraext_s = 1.0e08, & !: extinction coefficient of radiation in the snow 170 zkimin = 0.10 , & !: minimum ice thermal conductivity 171 zht_smin = 1.0e-4 !: minimum snow depth 172 173 REAL(wp) :: & ! local variables 174 ztmelt_i, & ! ice melting temperature 175 zerritmax ! current maximal error on temperature 176 177 REAL(wp), DIMENSION(jpij) :: & 178 zerrit, & ! current error on temperature 179 zdifcase, & ! case of the equation resolution (1->4) 180 zftrice, & ! solar radiation transmitted through the ice 181 zihic, zhsu 182 183 ! 184 !------------------------------------------------------------------------------! 185 ! 1) Initialization ! 186 !------------------------------------------------------------------------------! 187 ! 188 DO ji = kideb , kiut 189 ! is there snow or not 190 isnow(ji)= INT ( 1.0 - MAX( 0.0 , SIGN (1.0, - ht_s_b(ji) ) ) ) 191 ! surface temperature of fusion 192 ztfs(ji) = isnow(ji) * rtt + (1.0-isnow(ji)) * rtt 193 ! layer thickness 194 zh_i(ji) = ht_i_b(ji) / nlay_i 195 zh_s(ji) = ht_s_b(ji) / nlay_s 196 END DO 197 198 !-------------------- 199 ! Ice / snow layers 200 !-------------------- 201 202 z_s(:,0) = 0.0 ! vert. coord. of the up. lim. of the 1st snow layer 203 z_i(:,0) = 0.0 ! vert. coord. of the up. lim. of the 1st ice layer 204 205 DO layer = 1, nlay_s 206 DO ji = kideb , kiut 207 ! vert. coord of the up. lim. of the layer-th snow layer 208 z_s(ji,layer) = z_s(ji,layer-1) + ht_s_b(ji) / nlay_s 209 END DO 210 END DO 211 212 DO layer = 1, nlay_i 213 DO ji = kideb , kiut 214 ! vert. coord of the up. lim. of the layer-th ice layer 215 z_i(ji,layer) = z_i(ji,layer-1) + ht_i_b(ji) / nlay_i 216 END DO 217 END DO 218 ! 219 !------------------------------------------------------------------------------| 220 ! 2) Radiations | 221 !------------------------------------------------------------------------------| 222 ! 223 !------------------- 224 ! Computation of i0 225 !------------------- 226 ! i0 describes the fraction of solar radiation which does not contribute 227 ! to the surface energy budget but rather penetrates inside the ice. 228 ! We assume that no radiation is transmitted through the snow 229 ! If there is no no snow 230 ! zfsw = (1-i0).qsr_ice is absorbed at the surface 231 ! zftrice = io.qsr_ice is below the surface 232 ! fstbif = io.qsr_ice.exp(-k(h_i)) transmitted below the ice 233 234 DO ji = kideb , kiut 235 ! switches 236 isnow(ji) = INT ( 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_s_b(ji) ) ) ) 237 ! hs > 0, isnow = 1 238 zhsu(ji) = hnzst !threshold for the computation of i0 239 zihic(ji) = MAX( zzero , 1.0 - ( ht_i_b(ji) / zhsu(ji) ) ) 240 241 i0(ji) = ( 1.0 - isnow(ji) ) * & 242 ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) ) 243 !fr1_i0_1d = i0 for a thin ice surface 244 !fr1_i0_2d = i0 for a thick ice surface 245 ! a function of the cloud cover 246 ! 247 !i0(ji) = (1.0-FLOAT(isnow(ji)))*3.0/(100*ht_s_b(ji)+10.0) 248 !formula used in Cice 249 END DO 250 251 !------------------------------------------------------- 252 ! Solar radiation absorbed / transmitted at the surface 253 ! Derivative of the non solar flux 254 !------------------------------------------------------- 255 DO ji = kideb , kiut 256 257 ! Shortwave radiation absorbed at surface 258 zfsw(ji) = qsr_ice_1d(ji) * ( 1 - i0(ji) ) 259 260 ! Solar radiation transmitted below the surface layer 261 zftrice(ji)= qsr_ice_1d(ji) * i0(ji) 262 263 ! derivative of incoming nonsolar flux 264 dzf(ji) = dqns_ice_1d(ji) 265 266 END DO 267 268 !--------------------------------------------------------- 269 ! Transmission - absorption of solar radiation in the ice 270 !--------------------------------------------------------- 271 272 DO ji = kideb , kiut 273 ! Initialization 274 zradtr_s(ji,0) = zftrice(ji) ! radiation penetrating through snow 275 END DO 276 277 ! Radiation through snow 278 DO layer = 1, nlay_s 279 DO ji = kideb , kiut 280 ! radiation transmitted below the layer-th snow layer 281 zradtr_s(ji,layer) = zradtr_s(ji,0) * EXP ( - zraext_s * ( MAX ( 0.0 , & 282 z_s(ji,layer) ) ) ) 283 ! radiation absorbed by the layer-th snow layer 284 zradab_s(ji,layer) = zradtr_s(ji,layer-1) - zradtr_s(ji,layer) 285 END DO 286 END DO 287 288 ! Radiation through ice 289 DO ji = kideb , kiut 290 zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * isnow(ji) + & 291 zftrice(ji) * ( 1 - isnow(ji) ) 292 END DO 293 294 DO layer = 1, nlay_i 295 DO ji = kideb , kiut 296 ! radiation transmitted below the layer-th ice layer 297 zradtr_i(ji,layer) = zradtr_i(ji,0) * EXP ( - kappa_i * ( MAX ( 0.0 , & 298 z_i(ji,layer) ) ) ) 299 ! radiation absorbed by the layer-th ice layer 300 zradab_i(ji,layer) = zradtr_i(ji,layer-1) - zradtr_i(ji,layer) 301 END DO 302 END DO 303 304 ! Radiation transmitted below the ice 305 DO ji = kideb , kiut 306 fstbif_1d(ji) = fstbif_1d(ji) + & 307 zradtr_i(ji,nlay_i) * a_i_b(ji) / at_i_b(ji) 308 END DO 309 310 ! +++++ 311 ! just to check energy conservation 312 DO ji = kideb , kiut 313 zji = MOD( npb(ji) - 1, jpi ) + 1 314 zjj = ( npb(ji) - 1 ) / jpi + 1 315 fstroc(zji,zjj,jl) = & 316 zradtr_i(ji,nlay_i) 317 END DO 318 ! +++++ 319 320 DO layer = 1, nlay_i 321 DO ji = kideb , kiut 322 radab(ji,layer) = zradab_i(ji,layer) 323 END DO 324 END DO 325 326 327 ! 328 !------------------------------------------------------------------------------| 329 ! 3) Iterative procedure begins | 330 !------------------------------------------------------------------------------| 331 ! 332 ! Old surface temperature 333 DO ji = kideb, kiut 334 ztsuold(ji) = t_su_b(ji) ! temperature at the beg of iter pr. 335 ztsuoldit(ji) = t_su_b(ji) ! temperature at the previous iter 336 t_su_b(ji) = MIN(t_su_b(ji),ztfs(ji)-0.00001) !necessary 337 zerrit(ji) = 1000.0 ! initial value of error 338 END DO 339 !RB Min global ?? 340 341 ! Old snow temperature 342 DO layer = 1, nlay_s 343 DO ji = kideb , kiut 344 ztsold(ji,layer) = t_s_b(ji,layer) 345 END DO 346 END DO 347 348 ! Old ice temperature 349 DO layer = 1, nlay_i 350 DO ji = kideb , kiut 351 ztiold(ji,layer) = t_i_b(ji,layer) 352 END DO 353 END DO 354 355 nconv = 0 ! number of iterations 356 zerritmax = 1000.0 ! maximal value of error on all points 357 358 DO WHILE ((zerritmax > maxer_i_thd).AND.(nconv < nconv_i_thd)) 359 360 nconv = nconv+1 361 362 ! 363 !------------------------------------------------------------------------------| 364 ! 4) Sea ice thermal conductivity | 365 !------------------------------------------------------------------------------| 366 ! 367 IF ( thcon_i_swi .EQ. 0 ) THEN 368 ! Untersteiner (1964) formula 369 DO ji = kideb , kiut 370 ztcond_i(ji,0) = rcdic + zbeta*s_i_b(ji,1) / & 371 MIN(-zeps,t_i_b(ji,1)-rtt) 372 ztcond_i(ji,0) = MAX(ztcond_i(ji,0),zkimin) 373 END DO 374 ENDIF 375 376 IF ( thcon_i_swi .EQ. 1 ) THEN 377 ! Pringle et al formula included, 378 ! 2.11 + 0.09 S/T - 0.011.T 379 DO ji = kideb , kiut 380 ztcond_i(ji,0) = rcdic + 0.09*s_i_b(ji,1) / & 381 MIN(-zeps,t_i_b(ji,1)-rtt) - & 382 0.011* ( t_i_b(ji,1) - rtt ) 383 ztcond_i(ji,0) = MAX(ztcond_i(ji,0),zkimin) 384 END DO 385 ENDIF 386 387 IF ( thcon_i_swi .EQ. 0 ) THEN ! Untersteiner 388 DO layer = 1, nlay_i-1 389 DO ji = kideb , kiut 390 ztcond_i(ji,layer) = rcdic + zbeta*( s_i_b(ji,layer) & 391 + s_i_b(ji,layer+1) ) / MIN(-zeps, & 392 t_i_b(ji,layer)+t_i_b(ji,layer+1)-2.0*rtt) 393 ztcond_i(ji,layer) = MAX(ztcond_i(ji,layer),zkimin) 394 END DO 395 END DO 396 ENDIF 397 398 IF ( thcon_i_swi .EQ. 1 ) THEN ! Pringle 399 DO layer = 1, nlay_i-1 400 DO ji = kideb , kiut 401 ztcond_i(ji,layer) = rcdic + 0.09*( s_i_b(ji,layer) & 402 + s_i_b(ji,layer+1) ) / MIN(-zeps, & 403 t_i_b(ji,layer)+t_i_b(ji,layer+1)-2.0*rtt) - & 404 0.011* ( t_i_b(ji,layer) + t_i_b(ji,layer+1) - 2.0*rtt ) 405 ztcond_i(ji,layer) = MAX(ztcond_i(ji,layer),zkimin) 406 END DO 407 END DO 408 ENDIF 409 410 IF ( thcon_i_swi .EQ. 0 ) THEN ! Untersteiner 411 DO ji = kideb , kiut 412 ztcond_i(ji,nlay_i) = rcdic + zbeta*s_i_b(ji,nlay_i) / & 413 MIN(-zeps,t_bo_b(ji)-rtt) 414 ztcond_i(ji,nlay_i) = MAX(ztcond_i(ji,nlay_i),zkimin) 415 END DO 416 ENDIF 417 418 IF ( thcon_i_swi .EQ. 1 ) THEN ! Pringle 419 DO ji = kideb , kiut 420 ztcond_i(ji,nlay_i) = rcdic + 0.09*s_i_b(ji,nlay_i) / & 421 MIN(-zeps,t_bo_b(ji)-rtt) - & 422 0.011* ( t_bo_b(ji) - rtt ) 423 ztcond_i(ji,nlay_i) = MAX(ztcond_i(ji,nlay_i),zkimin) 424 END DO 425 ENDIF 426 ! 427 !------------------------------------------------------------------------------| 428 ! 5) kappa factors | 429 !------------------------------------------------------------------------------| 430 ! 431 DO ji = kideb, kiut 432 433 !-- Snow kappa factors 434 zkappa_s(ji,0) = rcdsn / MAX(zeps,zh_s(ji)) 435 zkappa_s(ji,nlay_s) = rcdsn / MAX(zeps,zh_s(ji)) 436 END DO 437 438 DO layer = 1, nlay_s-1 439 DO ji = kideb , kiut 440 zkappa_s(ji,layer) = 2.0 * rcdsn / & 441 MAX(zeps,2.0*zh_s(ji)) 442 END DO 443 END DO 444 445 DO layer = 1, nlay_i-1 446 DO ji = kideb , kiut 447 !-- Ice kappa factors 448 zkappa_i(ji,layer) = 2.0*ztcond_i(ji,layer)/ & 449 MAX(zeps,2.0*zh_i(ji)) 450 END DO 451 END DO 452 453 DO ji = kideb , kiut 454 zkappa_i(ji,0) = ztcond_i(ji,0)/MAX(zeps,zh_i(ji)) 455 zkappa_i(ji,nlay_i) = ztcond_i(ji,nlay_i) / MAX(zeps,zh_i(ji)) 456 !-- Interface 457 zkappa_s(ji,nlay_s) = 2.0*rcdsn*ztcond_i(ji,0)/MAX(zeps, & 458 (ztcond_i(ji,0)*zh_s(ji) + rcdsn*zh_i(ji))) 459 zkappa_i(ji,0) = zkappa_s(ji,nlay_s)*isnow(ji) & 460 + zkappa_i(ji,0)*(1.0-isnow(ji)) 461 END DO 462 ! 463 !------------------------------------------------------------------------------| 464 ! 6) Sea ice specific heat, eta factors | 465 !------------------------------------------------------------------------------| 466 ! 467 DO layer = 1, nlay_i 468 DO ji = kideb , kiut 469 ztitemp(ji,layer) = t_i_b(ji,layer) 470 zspeche_i(ji,layer) = cpic + zgamma*s_i_b(ji,layer)/ & 471 MAX((t_i_b(ji,layer)-rtt)*(ztiold(ji,layer)-rtt),zeps) 472 zeta_i(ji,layer) = rdt_ice / MAX(rhoic*zspeche_i(ji,layer)*zh_i(ji), & 473 zeps) 474 END DO 475 END DO 476 477 DO layer = 1, nlay_s 478 DO ji = kideb , kiut 479 ztstemp(ji,layer) = t_s_b(ji,layer) 480 zeta_s(ji,layer) = rdt_ice / MAX(rhosn*cpic*zh_s(ji),zeps) 481 END DO 482 END DO 483 ! 484 !------------------------------------------------------------------------------| 485 ! 7) surface flux computation | 486 !------------------------------------------------------------------------------| 487 ! 488 DO ji = kideb , kiut 489 490 ! update of the non solar flux according to the update in T_su 491 qnsr_ice_1d(ji) = qnsr_ice_1d(ji) + dqns_ice_1d(ji) * & 492 ( t_su_b(ji) - ztsuoldit(ji) ) 493 494 ! update incoming flux 495 zf(ji) = zfsw(ji) & ! net absorbed solar radiation 496 + qnsr_ice_1d(ji) ! non solar total flux 497 ! (LWup, LWdw, SH, LH) 498 499 END DO 500 501 ! 502 !------------------------------------------------------------------------------| 503 ! 8) tridiagonal system terms | 504 !------------------------------------------------------------------------------| 505 ! 506 !!layer denotes the number of the layer in the snow or in the ice 507 !!numeq denotes the reference number of the equation in the tridiagonal 508 !!system, terms of tridiagonal system are indexed as following : 509 !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one 510 511 !!ice interior terms (top equation has the same form as the others) 512 513 DO numeq=1,jkmax+2 514 DO ji = kideb , kiut 515 ztrid(ji,numeq,1) = 0. 516 ztrid(ji,numeq,2) = 0. 517 ztrid(ji,numeq,3) = 0. 518 zindterm(ji,numeq)= 0. 519 zindtbis(ji,numeq)= 0. 520 zdiagbis(ji,numeq)= 0. 521 ENDDO 522 ENDDO 523 524 DO numeq = nlay_s + 2, nlay_s + nlay_i 525 DO ji = kideb , kiut 526 layer = numeq - nlay_s - 1 527 ztrid(ji,numeq,1) = - zeta_i(ji,layer)*zkappa_i(ji,layer-1) 528 ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,layer)*(zkappa_i(ji,layer-1) + & 529 zkappa_i(ji,layer)) 530 ztrid(ji,numeq,3) = - zeta_i(ji,layer)*zkappa_i(ji,layer) 531 zindterm(ji,numeq) = ztiold(ji,layer) + zeta_i(ji,layer)* & 532 zradab_i(ji,layer) 533 END DO 534 ENDDO 535 536 numeq = nlay_s + nlay_i + 1 537 DO ji = kideb , kiut 543 538 !!ice bottom term 544 539 ztrid(ji,numeq,1) = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1) 545 540 ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,nlay_i)*( zkappa_i(ji,nlay_i)*zg1 & 546 541 + zkappa_i(ji,nlay_i-1) ) 547 542 ztrid(ji,numeq,3) = 0.0 548 543 zindterm(ji,numeq) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i)* & 549 550 551 552 553 554 544 ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i)*zg1 & 545 * t_bo_b(ji) ) 546 ENDDO 547 548 549 DO ji = kideb , kiut 555 550 IF ( ht_s_b(ji).gt.0.0 ) THEN 556 !557 !------------------------------------------------------------------------------|558 ! snow-covered cells |559 !------------------------------------------------------------------------------|560 !561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 !------------------------------------------------------------------------------|582 ! case 1 : no surface melting - snow present |583 !------------------------------------------------------------------------------|584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 !602 !------------------------------------------------------------------------------|603 ! case 2 : surface is melting - snow present |604 !------------------------------------------------------------------------------|605 !606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 !621 !------------------------------------------------------------------------------|622 ! cells without snow |623 !------------------------------------------------------------------------------|624 !625 626 !627 !------------------------------------------------------------------------------|628 ! case 3 : no surface melting - no snow |629 !------------------------------------------------------------------------------|630 !631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 ENDIF662 663 664 665 !666 !------------------------------------------------------------------------------|667 ! case 4 : surface is melting - no snow |668 !------------------------------------------------------------------------------|669 !670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 !699 !------------------------------------------------------------------------------|700 ! 9) tridiagonal system solving |701 !------------------------------------------------------------------------------|702 !703 704 ! Solve the tridiagonal system with Gauss elimination method.705 ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON,706 ! McGraw-Hill 1984.707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 ! ice temperatures730 731 732 733 734 735 736 737 738 739 740 741 551 ! 552 !------------------------------------------------------------------------------| 553 ! snow-covered cells | 554 !------------------------------------------------------------------------------| 555 ! 556 !!snow interior terms (bottom equation has the same form as the others) 557 DO numeq = 3, nlay_s + 1 558 layer = numeq - 1 559 ztrid(ji,numeq,1) = - zeta_s(ji,layer)*zkappa_s(ji,layer-1) 560 ztrid(ji,numeq,2) = 1.0 + zeta_s(ji,layer)*( zkappa_s(ji,layer-1) + & 561 zkappa_s(ji,layer) ) 562 ztrid(ji,numeq,3) = - zeta_s(ji,layer)*zkappa_s(ji,layer) 563 zindterm(ji,numeq) = ztsold(ji,layer) + zeta_s(ji,layer)* & 564 zradab_s(ji,layer) 565 END DO 566 567 !!case of only one layer in the ice (ice equation is altered) 568 IF ( nlay_i.eq.1 ) THEN 569 ztrid(ji,nlay_s+2,3) = 0.0 570 zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1)* & 571 t_bo_b(ji) 572 ENDIF 573 574 IF ( t_su_b(ji) .LT. rtt ) THEN 575 576 !------------------------------------------------------------------------------| 577 ! case 1 : no surface melting - snow present | 578 !------------------------------------------------------------------------------| 579 zdifcase(ji) = 1.0 580 numeqmin(ji) = 1 581 numeqmax(ji) = nlay_i + nlay_s + 1 582 583 !!surface equation 584 ztrid(ji,1,1) = 0.0 585 ztrid(ji,1,2) = dzf(ji) - zg1s*zkappa_s(ji,0) 586 ztrid(ji,1,3) = zg1s*zkappa_s(ji,0) 587 zindterm(ji,1) = dzf(ji)*t_su_b(ji) - zf(ji) 588 589 !!first layer of snow equation 590 ztrid(ji,2,1) = - zkappa_s(ji,0)*zg1s*zeta_s(ji,1) 591 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1)*(zkappa_s(ji,1) + zkappa_s(ji,0)*zg1s) 592 ztrid(ji,2,3) = - zeta_s(ji,1)* zkappa_s(ji,1) 593 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1)*zradab_s(ji,1) 594 595 ELSE 596 ! 597 !------------------------------------------------------------------------------| 598 ! case 2 : surface is melting - snow present | 599 !------------------------------------------------------------------------------| 600 ! 601 zdifcase(ji) = 2.0 602 numeqmin(ji) = 2 603 numeqmax(ji) = nlay_i + nlay_s + 1 604 605 !!first layer of snow equation 606 ztrid(ji,2,1) = 0.0 607 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + & 608 zkappa_s(ji,0) * zg1s ) 609 ztrid(ji,2,3) = - zeta_s(ji,1)*zkappa_s(ji,1) 610 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * & 611 ( zradab_s(ji,1) + & 612 zkappa_s(ji,0) * zg1s * t_su_b(ji) ) 613 ENDIF 614 ELSE 615 ! 616 !------------------------------------------------------------------------------| 617 ! cells without snow | 618 !------------------------------------------------------------------------------| 619 ! 620 IF (t_su_b(ji) .LT. rtt) THEN 621 ! 622 !------------------------------------------------------------------------------| 623 ! case 3 : no surface melting - no snow | 624 !------------------------------------------------------------------------------| 625 ! 626 zdifcase(ji) = 3.0 627 numeqmin(ji) = nlay_s + 1 628 numeqmax(ji) = nlay_i + nlay_s + 1 629 630 !!surface equation 631 ztrid(ji,numeqmin(ji),1) = 0.0 632 ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0)*zg1 633 ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0)*zg1 634 zindterm(ji,numeqmin(ji)) = dzf(ji)*t_su_b(ji) - zf(ji) 635 636 !!first layer of ice equation 637 ztrid(ji,numeqmin(ji)+1,1) = - zkappa_i(ji,0) * zg1 * zeta_i(ji,1) 638 ztrid(ji,numeqmin(ji)+1,2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) & 639 + zkappa_i(ji,0) * zg1 ) 640 ztrid(ji,numeqmin(ji)+1,3) = - zeta_i(ji,1)*zkappa_i(ji,1) 641 zindterm(ji,numeqmin(ji)+1)= ztiold(ji,1) + zeta_i(ji,1)*zradab_i(ji,1) 642 643 !!case of only one layer in the ice (surface & ice equations are altered) 644 645 IF (nlay_i.eq.1) THEN 646 ztrid(ji,numeqmin(ji),1) = 0.0 647 ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0)*2.0 648 ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0)*2.0 649 ztrid(ji,numeqmin(ji)+1,1) = -zkappa_i(ji,0)*2.0*zeta_i(ji,1) 650 ztrid(ji,numeqmin(ji)+1,2) = 1.0 + zeta_i(ji,1)*(zkappa_i(ji,0)*2.0 + & 651 zkappa_i(ji,1)) 652 ztrid(ji,numeqmin(ji)+1,3) = 0.0 653 654 zindterm(ji,numeqmin(ji)+1) = ztiold(ji,1) + zeta_i(ji,1)* & 655 ( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_b(ji) ) 656 ENDIF 657 658 ELSE 659 660 ! 661 !------------------------------------------------------------------------------| 662 ! case 4 : surface is melting - no snow | 663 !------------------------------------------------------------------------------| 664 ! 665 zdifcase(ji) = 4.0 666 numeqmin(ji) = nlay_s + 2 667 numeqmax(ji) = nlay_i + nlay_s + 1 668 669 !!first layer of ice equation 670 ztrid(ji,numeqmin(ji),1) = 0.0 671 ztrid(ji,numeqmin(ji),2) = 1.0 + zeta_i(ji,1)*(zkappa_i(ji,1) + zkappa_i(ji,0)* & 672 zg1) 673 ztrid(ji,numeqmin(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1) 674 zindterm(ji,numeqmin(ji)) = ztiold(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + & 675 zkappa_i(ji,0) * zg1 * t_su_b(ji) ) 676 677 !!case of only one layer in the ice (surface & ice equations are altered) 678 IF (nlay_i.eq.1) THEN 679 ztrid(ji,numeqmin(ji),1) = 0.0 680 ztrid(ji,numeqmin(ji),2) = 1.0 + zeta_i(ji,1)*(zkappa_i(ji,0)*2.0 + & 681 zkappa_i(ji,1)) 682 ztrid(ji,numeqmin(ji),3) = 0.0 683 zindterm(ji,numeqmin(ji)) = ztiold(ji,1) + zeta_i(ji,1)* & 684 (zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_b(ji)) & 685 + t_su_b(ji)*zeta_i(ji,1)*zkappa_i(ji,0)*2.0 686 ENDIF 687 688 ENDIF 689 ENDIF 690 691 END DO 692 693 ! 694 !------------------------------------------------------------------------------| 695 ! 9) tridiagonal system solving | 696 !------------------------------------------------------------------------------| 697 ! 698 699 ! Solve the tridiagonal system with Gauss elimination method. 700 ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, 701 ! McGraw-Hill 1984. 702 703 maxnumeqmax = 0 704 minnumeqmin = jkmax+4 705 706 DO ji = kideb , kiut 707 zindtbis(ji,numeqmin(ji)) = zindterm(ji,numeqmin(ji)) 708 zdiagbis(ji,numeqmin(ji)) = ztrid(ji,numeqmin(ji),2) 709 minnumeqmin = MIN(numeqmin(ji),minnumeqmin) 710 maxnumeqmax = MAX(numeqmax(ji),maxnumeqmax) 711 END DO 712 713 DO layer = minnumeqmin+1, maxnumeqmax 714 DO ji = kideb , kiut 715 numeq = min(max(numeqmin(ji)+1,layer),numeqmax(ji)) 716 zdiagbis(ji,numeq) = ztrid(ji,numeq,2) - ztrid(ji,numeq,1)* & 717 ztrid(ji,numeq-1,3)/zdiagbis(ji,numeq-1) 718 zindtbis(ji,numeq) = zindterm(ji,numeq) - ztrid(ji,numeq,1)* & 719 zindtbis(ji,numeq-1)/zdiagbis(ji,numeq-1) 720 END DO 721 END DO 722 723 DO ji = kideb , kiut 724 ! ice temperatures 725 t_i_b(ji,nlay_i) = zindtbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji)) 726 END DO 727 728 DO numeq = nlay_i + nlay_s + 1, nlay_s + 2, -1 729 DO ji = kideb , kiut 730 layer = numeq - nlay_s - 1 731 t_i_b(ji,layer) = (zindtbis(ji,numeq) - ztrid(ji,numeq,3)* & 732 t_i_b(ji,layer+1))/zdiagbis(ji,numeq) 733 END DO 734 END DO 735 736 DO ji = kideb , kiut 742 737 ! snow temperatures 743 738 IF (ht_s_b(ji).GT.0) & 744 t_s_b(ji,nlay_s) = (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) &745 746 739 t_s_b(ji,nlay_s) = (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) & 740 * t_i_b(ji,1))/zdiagbis(ji,nlay_s+1) & 741 * MAX(0.0,SIGN(1.0,ht_s_b(ji)-zeps)) 747 742 748 743 ! surface temperature … … 750 745 ztsuoldit(ji) = t_su_b(ji) 751 746 IF (t_su_b(ji) .LT. ztfs(ji)) & 752 t_su_b(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* &753 754 755 756 757 !758 !--------------------------------------------------------------------------759 ! 10) Has the scheme converged ?, end of the iterative procedure |760 !--------------------------------------------------------------------------761 !762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 747 t_su_b(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* & 748 ( isnow(ji)*t_s_b(ji,1) + & 749 (1.0-isnow(ji))*t_i_b(ji,1) ) ) / & 750 zdiagbis(ji,numeqmin(ji)) 751 END DO 752 ! 753 !-------------------------------------------------------------------------- 754 ! 10) Has the scheme converged ?, end of the iterative procedure | 755 !-------------------------------------------------------------------------- 756 ! 757 ! check that nowhere it has started to melt 758 ! zerrit(ji) is a measure of error, it has to be under maxer_i_thd 759 DO ji = kideb , kiut 760 t_su_b(ji) = MAX(MIN(t_su_b(ji),ztfs(ji)),190.0) 761 zerrit(ji) = ABS(t_su_b(ji)-ztsuoldit(ji)) 762 END DO 763 764 DO layer = 1, nlay_s 765 DO ji = kideb , kiut 766 zji = MOD( npb(ji) - 1, jpi ) + 1 767 zjj = ( npb(ji) - 1 ) / jpi + 1 768 t_s_b(ji,layer) = MAX(MIN(t_s_b(ji,layer),rtt),190.0) 769 zerrit(ji) = MAX(zerrit(ji),ABS(t_s_b(ji,layer) & 770 - ztstemp(ji,layer))) 771 END DO 772 END DO 773 774 DO layer = 1, nlay_i 775 DO ji = kideb , kiut 776 ztmelt_i = -tmut*s_i_b(ji,layer) +rtt 777 t_i_b(ji,layer) = MAX(MIN(t_i_b(ji,layer),ztmelt_i),190.0) 778 zerrit(ji) = MAX(zerrit(ji),ABS(t_i_b(ji,layer) - ztitemp(ji,layer))) 779 END DO 780 END DO 781 782 ! Compute spatial maximum over all errors 783 ! note that this could be optimized substantially by iterating only 784 ! the non-converging points 785 zerritmax = 0.0 786 DO ji = kideb , kiut 787 zerritmax = MAX(zerritmax,zerrit(ji)) 788 END DO 789 IF( lk_mpp ) CALL mpp_max(zerritmax, kcom=ncomm_ice) 795 790 796 791 END DO ! End of the do while iterative procedure … … 800 795 801 796 802 !803 !--------------------------------------------------------------------------804 ! 11) Fluxes at the interfaces |805 !--------------------------------------------------------------------------806 !807 808 ! update of latent heat fluxes809 810 811 812 ! surface ice conduction flux813 814 815 816 817 818 819 ! bottom ice conduction flux820 821 822 823 824 825 826 827 828 829 830 DO ji = kideb, kiut831 ! Upper snow value832 fc_s(ji,0) = - isnow(ji)* &833 834 835 ! Bott. snow value836 fc_s(ji,1) = - isnow(ji)* &837 838 839 END DO840 841 ! Upper ice layer842 DO ji = kideb, kiut843 fc_i(ji,0) = - isnow(ji) * & ! interface flux if there is snow844 845 846 847 END DO848 849 ! Internal ice layers850 DO layer = 1, nlay_i - 1851 DO ji = kideb, kiut852 fc_i(ji,layer) = - zkappa_i(ji,layer) * ( t_i_b(ji,layer+1) - &853 854 zji = MOD( npb(ji) - 1, jpi ) + 1855 zjj = ( npb(ji) - 1 ) / jpi + 1856 END DO857 END DO858 859 ! Bottom ice layers860 DO ji = kideb, kiut861 fc_i(ji,nlay_i) = - zkappa_i(ji,nlay_i)* &862 863 END DO864 865 866 867 797 ! 798 !-------------------------------------------------------------------------- 799 ! 11) Fluxes at the interfaces | 800 !-------------------------------------------------------------------------- 801 ! 802 DO ji = kideb, kiut 803 ! update of latent heat fluxes 804 qla_ice_1d (ji) = qla_ice_1d (ji) + & 805 dqla_ice_1d(ji) * ( t_su_b(ji) - ztsuold(ji) ) 806 807 ! surface ice conduction flux 808 isnow(ji) = int(1.0-max(0.0,sign(1.0,-ht_s_b(ji)))) 809 fc_su(ji) = - isnow(ji)*zkappa_s(ji,0)*zg1s*(t_s_b(ji,1) - & 810 t_su_b(ji)) & 811 - (1.0-isnow(ji))*zkappa_i(ji,0)*zg1* & 812 (t_i_b(ji,1) - t_su_b(ji)) 813 814 ! bottom ice conduction flux 815 fc_bo_i(ji) = - zkappa_i(ji,nlay_i)* & 816 ( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) ) 817 818 END DO 819 820 !-------------------------! 821 ! Heat conservation ! 822 !-------------------------! 823 IF ( con_i ) THEN 824 825 DO ji = kideb, kiut 826 ! Upper snow value 827 fc_s(ji,0) = - isnow(ji)* & 828 zkappa_s(ji,0) * zg1s * ( t_s_b(ji,1) - & 829 t_su_b(ji) ) 830 ! Bott. snow value 831 fc_s(ji,1) = - isnow(ji)* & 832 zkappa_s(ji,1) * ( t_i_b(ji,1) - & 833 t_s_b(ji,1) ) 834 END DO 835 836 ! Upper ice layer 837 DO ji = kideb, kiut 838 fc_i(ji,0) = - isnow(ji) * & ! interface flux if there is snow 839 ( zkappa_i(ji,0) * ( t_i_b(ji,1) - t_s_b(ji,nlay_s ) ) ) & 840 - ( 1.0 - isnow(ji) ) * ( zkappa_i(ji,0) * & 841 zg1 * ( t_i_b(ji,1) - t_su_b(ji) ) ) ! upper flux if not 842 END DO 843 844 ! Internal ice layers 845 DO layer = 1, nlay_i - 1 846 DO ji = kideb, kiut 847 fc_i(ji,layer) = - zkappa_i(ji,layer) * ( t_i_b(ji,layer+1) - & 848 t_i_b(ji,layer) ) 849 zji = MOD( npb(ji) - 1, jpi ) + 1 850 zjj = ( npb(ji) - 1 ) / jpi + 1 851 END DO 852 END DO 853 854 ! Bottom ice layers 855 DO ji = kideb, kiut 856 fc_i(ji,nlay_i) = - zkappa_i(ji,nlay_i)* & 857 ( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) ) 858 END DO 859 860 ENDIF 861 862 END SUBROUTINE lim_thd_dif 868 863 869 864 #else … … 876 871 END SUBROUTINE lim_thd_dif 877 872 #endif 878 873 END MODULE limthd_dif
Note: See TracChangeset
for help on using the changeset viewer.
