- Timestamp:
- 2017-09-18T16:54:04+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90
r8531 r8534 2 2 !!====================================================================== 3 3 !! *** MODULE icethd_zdf *** 4 !! heat diffusion in sea ice 5 !! computation of surface and inner T 4 !! sea-ice: vertical heat diffusion in sea ice (computation of temperatures) 6 5 !!====================================================================== 7 6 !! History : LIM ! 02-2003 (M. Vancoppenolle) original 1D code … … 14 13 #if defined key_lim3 15 14 !!---------------------------------------------------------------------- 16 !! 'key_lim3' LIM3sea-ice model15 !! 'key_lim3' ESIM sea-ice model 17 16 !!---------------------------------------------------------------------- 18 USE par_oce ! ocean parameters17 USE dom_oce ! ocean space and time domain 19 18 USE phycst ! physical constants (ocean directory) 20 19 USE ice ! sea-ice: variables 21 USE ice1D ! sea-ice: thermodynamics 20 USE ice1D ! sea-ice: thermodynamics variables 22 21 ! 23 22 USE in_out_manager ! I/O manager 24 23 USE lib_mpp ! MPP library 25 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)24 USE lib_fortran ! fortran utilities (glob_sum + no signed zero) 26 25 27 26 IMPLICIT NONE 28 27 PRIVATE 29 28 30 PUBLIC ice_thd_zdf ! called by ice _thd31 PUBLIC ice_thd_zdf_init ! called by ice _stp29 PUBLIC ice_thd_zdf ! called by icethd 30 PUBLIC ice_thd_zdf_init ! called by icestp 32 31 33 32 !!** namelist (namthd_zdf) ** … … 47 46 48 47 SUBROUTINE ice_thd_zdf 49 !!------------------------------------------------------------------ 48 !!------------------------------------------------------------------- 50 49 !! *** ROUTINE ice_thd_zdf *** 51 50 !! ** Purpose : … … 63 62 !! 64 63 !! The successive steps of this routine are 65 !! 1. Thermal conductivity at the interfaces of the ice layers 66 !! 2. Internal absorbed radiation 67 !! 3. Scale factors due to non-uniform grid 64 !! 1. initialization of ice-snow layers thicknesses 65 !! 2. Internal absorbed and transmitted radiation 66 !! Then iterative procedure begins 67 !! 3. Thermal conductivity 68 68 !! 4. Kappa factors 69 !! Then iterative procedure begins70 69 !! 5. specific heat in the ice 71 70 !! 6. eta factors … … 75 74 !! Iterative procedure ends according to a criterion on evolution 76 75 !! of temperature 76 !! 10. Fluxes at the interfaces 77 77 !! 78 78 !! ** Inputs / Ouputs : (global commons) … … 82 82 !! number of layers in the ice/snow: nlay_i, nlay_s 83 83 !! total ice/snow thickness : ht_i_1d, ht_s_1d 84 !!------------------------------------------------------------------ 84 !!------------------------------------------------------------------- 85 85 INTEGER :: ji, jk ! spatial loop index 86 86 INTEGER :: numeq ! current reference number of equation … … 149 149 END DO 150 150 151 !------------------ ------------------------------------------------------------!152 ! 1) Initialization !153 !------------------ ------------------------------------------------------------!151 !------------------ 152 ! 1) Initialization 153 !------------------ 154 154 DO ji = 1, nidx 155 155 isnow(ji)= 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) ) ! is there snow or not … … 173 173 t_su_1d(1:nidx) = MIN( t_su_1d(1:nidx), rt0 - ztsu_err ) ! necessary 174 174 ! 175 !------------- -----------------------------------------------------------------|176 ! 2) Radiation |177 !------------- -----------------------------------------------------------------|175 !------------- 176 ! 2) Radiation 177 !------------- 178 178 ! 179 179 z1_hsu = 1._wp / 0.1_wp ! threshold for the computation of i0 180 180 DO ji = 1, nidx 181 !------------------- 182 ! Computation of i0 183 !------------------- 181 ! --- Computation of i0 --- ! 184 182 ! i0 describes the fraction of solar radiation which does not contribute 185 183 ! to the surface energy budget but rather penetrates inside the ice. … … 193 191 i0(ji) = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zfac * fr2_i0_1d(ji) ) 194 192 195 !------------------------------------------------------- 196 ! Solar radiation absorbed / transmitted at the surface 197 ! Derivative of the non solar flux 198 !------------------------------------------------------- 193 ! --- Solar radiation absorbed / transmitted at the surface --- ! 194 ! Derivative of the non solar flux 199 195 zfsw (ji) = qsr_ice_1d(ji) * ( 1 - i0(ji) ) ! Shortwave radiation absorbed at surface 200 196 zftrice(ji) = qsr_ice_1d(ji) * i0(ji) ! Solar radiation transmitted below the surface layer … … 203 199 END DO 204 200 205 !--------------------------------------------------------- 206 ! Transmission - absorption of solar radiation in the ice 207 !--------------------------------------------------------- 201 ! --- Transmission/absorption of solar radiation in the ice --- ! 208 202 zradtr_s(1:nidx,0) = zftrice(1:nidx) 209 203 DO jk = 1, nlay_s … … 227 221 228 222 ftr_ice_1d(1:nidx) = zradtr_i(1:nidx,nlay_i) ! record radiation transmitted below the ice 229 230 !------------------------------------------------------------------------------|231 ! 3) Iterative procedure begins |232 !------------------------------------------------------------------------------|233 223 ! 234 224 iconv = 0 ! number of iterations 235 225 zdti_max = 1000._wp ! maximal value of error on all points 236 DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max ) 226 ! !----------------------------! 227 DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max ) ! Iterative procedure begins ! 228 ! !----------------------------! 237 229 ! 238 230 iconv = iconv + 1 … … 241 233 ztsb(1:nidx,:) = t_s_1d(1:nidx,:) 242 234 ! 243 !------------------------------------------------------------------------------| 244 ! 4) Sea ice thermal conductivity | 245 !------------------------------------------------------------------------------| 246 ! 235 !-------------------------------- 236 ! 3) Sea ice thermal conductivity 237 !-------------------------------- 247 238 IF( ln_cndi_U64 ) THEN !-- Untersteiner (1964) formula: k = k0 + beta.S/T 248 239 ! … … 277 268 ztcond_i(1:nidx,:) = MAX( zkimin, ztcond_i(1:nidx,:) ) 278 269 ! 279 !------------------------------------------------------------------------------| 280 ! 5) G(he) - enhancement of thermal conductivity in mono-category case | 281 !------------------------------------------------------------------------------| 282 ! 270 !--- G(he) : enhancement of thermal conductivity in mono-category case 283 271 ! Computation of effective thermal conductivity G(h) 284 272 ! Used in mono-category case only to simulate an ITD implicitly … … 309 297 END SELECT 310 298 ! 311 !------------------------------------------------------------------------------| 312 ! 6) kappa factors | 313 !------------------------------------------------------------------------------| 314 ! 299 !----------------- 300 ! 4) kappa factors 301 !----------------- 315 302 !--- Snow 316 303 DO jk = 0, nlay_s-1 … … 338 325 END DO 339 326 ! 340 !------------------------------------------------------------------------------| 341 ! 7) Sea ice specific heat, eta factors | 342 !------------------------------------------------------------------------------| 343 ! 327 !-------------------------------------- 328 ! 5) Sea ice specific heat, eta factors 329 !-------------------------------------- 344 330 DO jk = 1, nlay_i 345 331 DO ji = 1, nidx … … 355 341 END DO 356 342 ! 357 !------------------------------------------------------------------------------| 358 ! 8) surface flux computation | 359 !------------------------------------------------------------------------------| 360 ! 343 !---------------------------- 344 ! 6) surface flux computation 345 !---------------------------- 361 346 IF ( ln_dqns_i ) THEN 362 347 DO ji = 1, nidx … … 370 355 END DO 371 356 ! 372 !------------------------------------------------------------------------------| 373 ! 9) tridiagonal system terms | 374 !------------------------------------------------------------------------------| 375 ! 357 !---------------------------- 358 ! 7) tridiagonal system terms 359 !---------------------------- 376 360 !!layer denotes the number of the layer in the snow or in the ice 377 361 !!numeq denotes the reference number of the equation in the tridiagonal … … 414 398 415 399 DO ji = 1, nidx 416 IF ( ht_s_1d(ji) > 0.0 ) THEN 417 ! 418 !------------------------------------------------------------------------------| 419 ! snow-covered cells | 420 !------------------------------------------------------------------------------| 400 ! !---------------------! 401 IF ( ht_s_1d(ji) > 0.0 ) THEN ! snow-covered cells ! 402 ! !---------------------! 421 403 ! 422 404 !!snow interior terms (bottom equation has the same form as the others) … … 435 417 ENDIF 436 418 437 IF ( t_su_1d(ji) < rt0 ) THEN 438 439 !------------------------------------------------------------------------------| 440 ! case 1 : no surface melting - snow present | 441 !------------------------------------------------------------------------------| 419 IF ( t_su_1d(ji) < rt0 ) THEN !-- case 1 : no surface melting 420 442 421 numeqmin(ji) = 1 443 422 numeqmax(ji) = nlay_i + nlay_s + 1 … … 455 434 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1) 456 435 457 ELSE 458 ! 459 !------------------------------------------------------------------------------| 460 ! case 2 : surface is melting - snow present | 461 !------------------------------------------------------------------------------| 436 ELSE !-- case 2 : surface is melting 462 437 ! 463 438 numeqmin(ji) = 2 … … 471 446 & ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 472 447 ENDIF 473 ELSE 448 ! !---------------------! 449 ELSE ! cells without snow ! 450 ! !---------------------! 474 451 ! 475 !------------------------------------------------------------------------------| 476 ! cells without snow | 477 !------------------------------------------------------------------------------| 478 ! 479 IF ( t_su_1d(ji) < rt0 ) THEN 480 ! 481 !------------------------------------------------------------------------------| 482 ! case 3 : no surface melting - no snow | 483 !------------------------------------------------------------------------------| 452 IF ( t_su_1d(ji) < rt0 ) THEN !-- case 1 : no surface melting 484 453 ! 485 454 numeqmin(ji) = nlay_s + 1 … … 512 481 ENDIF 513 482 514 ELSE 515 516 ! 517 !------------------------------------------------------------------------------| 518 ! case 4 : surface is melting - no snow | 519 !------------------------------------------------------------------------------| 520 ! 483 ELSE !-- case 2 : surface is melting 484 521 485 numeqmin(ji) = nlay_s + 2 522 486 numeqmax(ji) = nlay_i + nlay_s + 1 … … 543 507 END DO 544 508 ! 545 !------------------------------------------------------------------------------| 546 ! 10) tridiagonal system solving | 547 !------------------------------------------------------------------------------| 548 ! 509 !------------------------------ 510 ! 8) tridiagonal system solving 511 !------------------------------ 549 512 ! Solve the tridiagonal system with Gauss elimination method. 550 ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, 551 ! McGraw-Hill 1984. 513 ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984. 552 514 553 515 maxnumeqmax = 0 … … 595 557 END DO 596 558 ! 597 !-------------------------------------------------------------------------- 598 ! 11) Has the scheme converged ?, end of the iterative procedure | 599 !-------------------------------------------------------------------------- 600 ! 559 !-------------------------------------------------------------- 560 ! 9) Has the scheme converged ?, end of the iterative procedure 561 !-------------------------------------------------------------- 601 562 ! check that nowhere it has started to melt 602 563 ! zdti_max is a measure of error, it has to be under zdti_bnd … … 632 593 WRITE(numout,*) ' iconv : ', iconv 633 594 ENDIF 634 635 ! 636 !-------------------------------------------------------------------------! 637 ! 12) Fluxes at the interfaces ! 638 !-------------------------------------------------------------------------! 595 ! 596 !----------------------------- 597 ! 10) Fluxes at the interfaces 598 !----------------------------- 639 599 DO ji = 1, nidx 640 600 ! ! surface ice conduction flux … … 694 654 695 655 SUBROUTINE ice_thd_enmelt 696 !!------------------------------------------------------------------- ----656 !!------------------------------------------------------------------- 697 657 !! *** ROUTINE ice_thd_enmelt *** 698 658 !! … … 764 724 ! 765 725 IF ( ( ln_cndi_U64 .AND. ln_cndi_P07 ) .OR. ( .NOT.ln_cndi_U64 .AND. .NOT.ln_cndi_P07 ) ) THEN 766 CALL ctl_stop( 'ice_thd_zdf_init: choose one and only one formulation for thermal conducti vity(ln_cndi_U64 or ln_cndi_P07)' )726 CALL ctl_stop( 'ice_thd_zdf_init: choose one and only one formulation for thermal conduction (ln_cndi_U64 or ln_cndi_P07)' ) 767 727 ENDIF 768 728 ! … … 771 731 #else 772 732 !!---------------------------------------------------------------------- 773 !! Dummy ModuleNo ESIM sea-ice model733 !! Default option Dummy Module No ESIM sea-ice model 774 734 !!---------------------------------------------------------------------- 775 735 #endif
Note: See TracChangeset
for help on using the changeset viewer.
