Changeset 834 for trunk/NEMO/LIM_SRC_3/limthd_dh.F90
- Timestamp:
- 2008-03-07T18:11:35+01:00 (16 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/LIM_SRC_3/limthd_dh.F90 (modified) (43 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/LIM_SRC_3/limthd_dh.F90
r825 r834 1 1 MODULE limthd_dh 2 2 #if defined key_lim3 3 !!---------------------------------------------------------------------- 4 !! 'key_lim3' LIM3 sea-ice model 5 !!---------------------------------------------------------------------- 3 6 !!====================================================================== 4 7 !! *** MODULE limthd_dh *** … … 19 22 USE ice 20 23 USE par_ice 21 USE limicepoints22 24 23 25 IMPLICIT NONE … … 36 38 37 39 !!---------------------------------------------------------------------- 38 !! LIM 3.0, UCL-ASTR-LOCEAN-IPSL (200 5)40 !! LIM 3.0, UCL-ASTR-LOCEAN-IPSL (2008) 39 41 !!---------------------------------------------------------------------- 40 42 … … 70 72 !! ** References : Bitz and Lipscomb, JGR 99 71 73 !! Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646 72 !! Vancoppenolle Fichefet and Bitz, GRL 2005 74 !! Vancoppenolle, Fichefet and Bitz, GRL 2005 75 !! Vancoppenolle et al., OM08 73 76 !! 74 77 !! ** History : 75 78 !! original code 01-04 (LIM) 76 79 !! original routine 77 !! (05-2003) Martin Vancoppenolle, Louvain-La-Neuve, Belgium 78 !! (06/07-2005) Martin Vancoppenolle for the 3D version 80 !! (05-2003) M. Vancoppenolle, Louvain-La-Neuve, Belgium 81 !! (06/07-2005) 3D version 82 !! (03-2008) Clean code 79 83 !! 80 84 !!------------------------------------------------------------------ … … 90 94 jk , & !: ice layer index 91 95 isnow , & !: switch for presence (1) or absence (0) of snow 92 index , & !: ice point index (limicepoints)93 96 zji , & !: 2D corresponding indices to ji 94 97 zjj , & … … 107 110 zhnfi , & 108 111 zihg , & 109 zdhgm , &110 112 zdhnm , & 111 113 zhnnew , & 112 zdfseqv , &113 114 zeps = 1.0e-13, & 114 115 zhisn , & … … 116 117 !: entrapment 117 118 zds , & !: increment of bottom ice salinity 118 zoldsinew , & !: dummy argument for error diagnosis119 119 zcoeff , & !: dummy argument for snowfall partitioning 120 120 !: over ice and leads 121 zjiindex , & !: zdhs_temp, &122 121 zsm_snowice, & !: snow-ice salinity 123 122 zswi1 , & !: switch for computation of bottom salinity … … 159 158 160 159 REAL(wp), DIMENSION(jpij,jkmax) :: & 161 zqt_i_lay , & !: total ice heat content 162 zqt_s_lay !: total snow heat content 160 zqt_i_lay !: total ice heat content 163 161 164 162 ! Heat conservation … … 181 179 WRITE(numout,*) 'lim_thd_dh : computation of vertical snow/ice accretion/ablation' 182 180 WRITE(numout,*) '~~~~~~~~~' 181 183 182 zfsalt_melt(:) = 0.0 184 183 ftotal_fin(:) = 0.0 … … 211 210 dsm_i_se_1d(:) = 0.0 212 211 dsm_i_si_1d(:) = 0.0 213 ! ! heat conservation 214 ! IF ( con_i ) THEN 215 ! IF (jiindex_1d .GT.0) WRITE(numout,*) ' zf_surf : ', z_f_surf(jiindex_1d) 216 ! IF (jiindex_1d .GT.0) WRITE(numout,*) ' qnsr_ice : ', qnsr_ice_1d(jiindex_1d) 217 ! IF (jiindex_1d .GT.0) WRITE(numout,*) ' qsr_ice : ', qsr_ice_1d(jiindex_1d) 218 ! IF (jiindex_1d .GT.0) WRITE(numout,*) ' fc_su : ', fc_su(jiindex_1d) 219 ! IF (jiindex_1d .GT.0) WRITE(numout,*) ' i0 : ', i0(jiindex_1d) 220 ! ENDIF 221 212 ! 222 213 !------------------------------------------------------------------------------! 223 214 ! 2) Computing layer thicknesses and snow and sea-ice enthalpies. ! 224 215 !------------------------------------------------------------------------------! 225 226 !----------------- 216 ! 227 217 ! Layer thickness 228 !-----------------229 218 DO ji = kideb,kiut 230 219 zh_i(ji) = ht_i_b(ji) / nlay_i … … 232 221 END DO 233 222 234 !-------------------------------------- 235 ! Snow energy of melting, heat content 236 !-------------------------------------- 223 ! Total enthalpy of the snow 237 224 zqt_s(:) = 0.0 238 225 DO jk = 1, nlay_s 239 226 DO ji = kideb,kiut 240 ! this line could be removed241 ! q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus )242 ! heat conservation243 227 zqt_s(ji) = zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s 244 228 END DO 245 229 END DO 246 230 247 ! ! heat conservation 248 ! qt_s_in(:,jl) = zqt_s(:) 249 ! IF (jiindex_1d.GT.0) & 250 ! WRITE(numout,*) 'jl : ', jl, ' zqts : ', zqt_s(jiindex_1d) / & 251 ! rdt_ice 252 253 ! IF (jiindex_1d.GT.0) THEN 254 ! WRITE(numout,*) ' Vertical profile : ' 255 ! WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d) 256 ! WRITE(numout,*) ' qm0 : ', q_s_b(jiindex_1d,1:nlay_s) * ht_s_b(jiindex_1d) / & 257 ! rdt_ice 258 ! WRITE(numout,*) ' q_s_b : ', q_s_b(jiindex_1d,1) 259 ! WRITE(numout,*) ' t_s_b : ', t_s_b(jiindex_1d,1) 260 ! WRITE(numout,*) ' zthick0: ', ht_s_b(jiindex_1d) 261 ! ENDIF 262 263 !------------------------------------- 264 ! Ice energy of melting, heat content 265 !------------------------------------- 231 ! Total enthalpy of the ice 266 232 zqt_i(:) = 0.0 267 233 DO jk = 1, nlay_i 268 234 DO ji = kideb,kiut 269 ! Melting point in K270 ! ztmelts = - tmut * s_i_b(ji,jk) + rtt271 ! this line could be removed272 ! q_i_b(ji,jk) = rhoic * &273 ! ( cpic * ( ztmelts - t_i_b(ji,jk) ) &274 ! + lfus * ( 1.0 - (ztmelts-rtt) / &275 ! MIN( ( t_i_b(ji,jk) - rtt ) , - zeps ) ) &276 ! - rcp * ( ztmelts-rtt ) )277 235 zqt_i(ji) = zqt_i(ji) + q_i_b(ji,jk) * ht_i_b(ji) / nlay_i 278 236 zqt_i_lay(ji,jk) = q_i_b(ji,jk) * ht_i_b(ji) / nlay_i 279 237 END DO 280 238 END DO 281 ! IF (jiindex_1d.GT.0) &282 ! WRITE(numout,*) 'jl : ', jl, ' zqti : ', zqt_i(jiindex_1d) / &283 ! rdt_ice284 285 ! IF (jiindex_1d.GT.0) THEN286 ! WRITE(numout,*) ' --- Vertical profile : '287 ! WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d)288 ! WRITE(numout,*) ' ht_i_b : ', ht_i_b(jiindex_1d)289 ! WRITE(numout,*) ' q_i_b : ', q_i_b(jiindex_1d, 1:nlay_i)290 ! WRITE(numout,*) ' t_i_b : ', t_i_b(jiindex_1d, 1:nlay_i)291 ! WRITE(numout,*) ' qm0 : ', zqt_i_lay(jiindex_1d,1:nlay_i) / rdt_ice292 ! WRITE(numout,*) ' zthick0: ', ht_i_b(jiindex_1d) / nlay_i293 ! ENDIF294 295 239 ! 296 240 !------------------------------------------------------------------------------| … … 298 242 !------------------------------------------------------------------------------| 299 243 ! 300 !--------------------- 301 ! Snow precips / melt302 !--------------------- 244 !------------------------- 245 ! 3.1 Snow precips / melt 246 !------------------------- 303 247 ! Snow accumulation in one thermodynamic time step 304 248 ! snowfall is partitionned between leads and ice … … 306 250 ! but because of the winds, more snow falls on leads than on sea ice 307 251 ! and a greater fraction (1-at_i)^beta of the total mass of snow 308 ! (beta < 1) falls in leads 252 ! (beta < 1) falls in leads. 309 253 ! In reality, beta depends on wind speed, 310 254 ! and should decrease with increasing wind speed but here, it is 311 ! considered as a constant 312 ! an average value is 0.66 255 ! considered as a constant. an average value is 0.66 313 256 ! Martin Vancoppenolle, December 2006 314 257 … … 322 265 ! Melt of fallen snow 323 266 DO ji = kideb, kiut 324 325 267 ! tatm_ice is now in K 326 268 zqprec(ji) = rhosn * ( cpic * ( rtt - tatm_ice_1d(ji) ) + lfus ) … … 334 276 qt_s_in(ji,jl) = qt_s_in(ji,jl) + zqprec(ji) * zdh_s_pre(ji) 335 277 zqt_s(ji) = zqt_s(ji) + zqprec(ji) * zdh_s_pre(ji) 336 337 END DO 338 339 ! IF (jiindex_1d .GT. 0) THEN 340 ! WRITE(numout,*) ' zqfont_su : ', zqfont_su(jiindex_1d) 341 ! WRITE(numout,*) ' zqprec : ', zqprec(jiindex_1d)*( zdh_s_pre(jiindex_1d) + zdh_s_mel(jiindex_1d) ) / rdt_ice 342 ! WRITE(numout,*) ' tatm : ', tatm_ice_1d(jiindex_1d) 343 ! WRITE(numout,*) 'jl : ', jl, ' zqts : ', zqt_s(jiindex_1d) / & 344 ! rdt_ice 345 ! ENDIF 346 347 ! updated total snow heat content 278 END DO 279 280 ! Update total snow heat content 348 281 zqt_s(ji) = MAX ( zqt_s(ji) - zqfont_su(ji) , 0.0 ) 349 350 ! IF (jiindex_1D .GT. 0) &351 ! WRITE(numout,*) 'jl : ', jl, ' zqts : ', zqt_s(jiindex_1d) / &352 ! rdt_ice353 282 354 283 ! Snow melt due to surface heat imbalance … … 358 287 zqfont_su(ji) = MAX( 0.0 , - zh_s(ji) - zdeltah(ji,jk) ) * & 359 288 q_s_b(ji,jk) 360 ! IF (ji.eq.jiindex_1d) WRITE(numout,*) ' zqfont_su : ', &361 ! zqfont_su(jiindex_1d)362 289 zdeltah(ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) ) 363 290 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) !resulting melt of snow 364 291 END DO 365 292 END DO 366 367 ! !+++++368 ! IF (jiindex_1d .GT. 0 ) THEN369 ! WRITE(numout,*) ' dhs_pre :', zdh_s_pre(jiindex_1d)370 ! WRITE(numout,*) ' dhs_mel :', zdh_s_mel(jiindex_1d)371 ! WRITE(numout,*) ' zqfont_su : ', zqfont_su(jiindex_1d)372 ! ENDIF373 ! !+++++374 293 375 294 ! Apply snow melt to snow depth … … 389 308 END DO ! ji 390 309 391 ! !+++++ 392 ! IF (jiindex_1d .GT. 1) WRITE(numout,*) ' dhs_tot :', dh_s_tot(jiindex_1d) 393 ! IF (jiindex_1d .GT. 1) WRITE(numout,*) ' zqfont_su :', zqfont_su(jiindex_1d) 394 ! !+++++ 395 396 !---------------------- 397 ! Surface ice ablation 398 !---------------------- 310 !-------------------------- 311 ! 3.2 Surface ice ablation 312 !-------------------------- 399 313 DO ji = kideb, kiut 400 ! IF (ji.eq.jiindex_1d) WRITE(numout,*) ' zqfont_su : ', &401 ! zqfont_su(jiindex_1d)402 314 dh_i_surf(ji) = 0.0 403 ! Heat conservation test315 ! For heat conservation test 404 316 z_f_surf(ji) = zqfont_su(ji) / rdt_ice ! heat conservation test 405 317 zdq_i(ji) = 0.0 … … 408 320 DO jk = 1, nlay_i 409 321 DO ji = kideb, kiut 322 ! melt of layer jk 410 323 zdeltah(ji,jk) = - zqfont_su(ji) / q_i_b(ji,jk) 324 ! recompute heat available 411 325 zqfont_su(ji) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * & 412 326 q_i_b(ji,jk) 327 ! melt of layer jk cannot be higher than its thickness 413 328 zdeltah(ji,jk) = MAX( zdeltah(ji,jk) , - zh_i(ji) ) 329 ! update surface melt 414 330 dh_i_surf(ji) = dh_i_surf(ji) + zdeltah(ji,jk) 331 ! for energy conservation 415 332 zdq_i(ji) = zdq_i(ji) + zdeltah(ji,jk) * & 416 333 q_i_b(ji,jk) / rdt_ice 417 ! Contribution to salt flux -> to change334 ! contribution to ice-ocean salt flux 418 335 zji = MOD( npb(ji) - 1, jpi ) + 1 419 336 zjj = ( npb(ji) - 1 ) / jpi + 1 420 337 zfsalt_melt(ji) = zfsalt_melt(ji) + & 421 ! ( sss_io(zji,zjj) - s_i_b(ji,jk) ) * &422 338 ( sss_io(zji,zjj) - sm_i_b(ji) ) * & 423 339 a_i_b(ji) * & 424 340 MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice 425 !+++++++++++426 ! IF ( (zji.eq.jiindex) .AND. (zjj.eq.jjindex) ) THEN427 ! WRITE(numout,*) ' surf abl '428 ! WRITE(numout,*) ' zfsalt_melt : ', zfsalt_melt(ji) / ( sss_io(zji,zjj)+epsi16)429 ! WRITE(numout,*) ' sss_io : ', sss_io(zji,zjj)430 ! ENDIF431 !+++++++++++432 341 END DO ! ji 433 342 END DO ! jk … … 446 355 ENDIF 447 356 448 IF ( z_f_surf(ji) + zdq_i(ji) .GE. 1.0e-3 ) THEN!.GE. 1.0e-3 ) .AND. & 449 ! ( ( ht_i_b(ji) + dh_i_bott(ji) ) .GE. 1.0e-6 ) ) THEN 357 IF ( z_f_surf(ji) + zdq_i(ji) .GE. 1.0e-3 ) THEN! 450 358 WRITE(numout,*) ' ALERTE heat loss for surface melt ' 451 359 WRITE(numout,*) ' zji, zjj, jl :', zji, zjj, jl … … 471 379 ENDIF ! con_i 472 380 473 !------------------ 474 ! Snow sublimation 475 !------------------ 476 ! !++++++ 477 ! zqt_dummy(:) = 0.0 478 ! DO jk = 1, nlay_s 479 ! DO ji = kideb,kiut 480 ! q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) 481 ! ! heat conservation 482 ! zqt_dummy(ji) = zqt_dummy(ji) + q_s_b(ji,jk) * old_ht_s_b(ji) / nlay_s 483 ! END DO 484 ! END DO 485 486 ! IF (jiindex_1d .GT. 0) & 487 ! WRITE(numout,*) 'jl : ', jl, ' Old zqts : ', zqt_dummy(jiindex_1d) / & 488 ! rdt_ice 489 !++++++ 490 ! !++++++ 491 ! zqt_dummy(:) = 0.0 492 ! DO jk = 1, nlay_s 493 ! DO ji = kideb,kiut 494 ! q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) 495 ! ! heat conservation 496 ! zqt_dummy(ji) = zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s 497 ! END DO 498 ! END DO 499 ! IF (jiindex_1d .GT. 0) & 500 ! WRITE(numout,*) 'jl : ', jl, ' New zqts : ', zqt_dummy(jiindex_1d) / & 501 ! rdt_ice 502 ! !++++++ 381 !---------------------- 382 ! 3.3 Snow sublimation 383 !---------------------- 503 384 504 385 DO ji = kideb, kiut … … 510 391 ht_s_b(ji) = MAX( zzero , zdhcf ) 511 392 ! we recompute dh_s_tot 512 ! which may be different in case of total snow melt513 393 dh_s_tot(ji) = ht_s_b(ji) - zhsold(ji) 514 ! IF ( ht_s_b(ji) .LE. 0.0 ) THEN515 ! dh_s_tot(ji) = MAX( 0.0 , dh_s_tot(ji) )516 ! ENDIF517 518 394 qt_s_in(ji,jl) = qt_s_in(ji,jl) + zdh_s_sub(ji)*q_s_b(ji,1) 519 395 END DO !ji 520 396 521 ! IF ( jiindex_1d .GT. 0 ) THEN522 ! WRITE(numout,*) ' dh_s_sub : ', zdh_s_sub(jiindex_1d)523 ! WRITE(numout,*) ' dhs_tot : ', dh_s_tot(jiindex_1d)524 ! ENDIF525 526 !++++++527 397 zqt_dummy(:) = 0.0 528 398 DO jk = 1, nlay_s … … 533 403 END DO 534 404 END DO 535 ! IF (jiindex_1d .GT. 0) &536 ! WRITE(numout,*) 'jl : ', jl, ' Old zqts : ', zqt_dummy(jiindex_1d) / &537 ! rdt_ice538 ! !++++++539 540 ! IF (jiindex_1d .GT. 0) THEN541 ! WRITE(numout,*) ' t_s_b : ', t_s_b(jiindex_1d,1)542 ! WRITE(numout,*) ' q_s_b : ', q_s_b(jiindex_1d,1)543 ! ENDIF544 405 545 406 DO jk = 1, nlay_s !n … … 553 414 END DO 554 415 555 ! IF (jiindex_1d .GT. 0) THEN556 ! WRITE(numout,*) ' t_s_b : ', t_s_b(jiindex_1d,1)557 ! WRITE(numout,*) ' q_s_b : ', q_s_b(jiindex_1d,1)558 ! ENDIF559 560 416 ! 561 417 !------------------------------------------------------------------------------! … … 567 423 ! the inner conductive flux (fc_bo_i), from the open water heat flux 568 424 ! (qlbbqb) and the turbulent ocean flux (fbif). 569 ! fc_bo_i is positive downwards 570 ! fbif and qlbbq are positive to the ice 571 572 !--------------------------------------------- 573 ! Basal growth - salinity not varying in time 574 !--------------------------------------------- 575 IF ( num_sal .NE. 2 ) THEN 425 ! fc_bo_i is positive downwards. fbif and qlbbq are positive to the ice 426 427 !----------------------------------------------------- 428 ! 4.1 Basal growth - (a) salinity not varying in time 429 !----------------------------------------------------- 430 IF ( ( num_sal .NE. 2 ) .AND. ( num_sal .NE. 4 ) ) THEN 576 431 DO ji = kideb, kiut 577 432 IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0 ) THEN … … 590 445 dh_i_bott(ji) = - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + & 591 446 qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 592 ! IF ( ji .EQ. jiindex_1d ) THEN593 ! WRITE(numout,*) ' s_i_new : ', s_i_new(ji)594 ! WRITE(numout,*) ' ztmelts : ', ztmelts595 ! WRITE(numout,*) ' t_bo : ', t_bo_b(ji)596 ! WRITE(numout,*) ' q_i_b : ', q_i_b(ji,nlay_i+1)597 ! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)598 ! ENDIF599 447 ENDIF ! heat budget 600 448 END DO ! ji 601 449 ENDIF ! num_sal 602 450 603 ! IF ( jiindex_1d .GT. 0 ) THEN 604 ! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(jiindex_1d) 605 ! WRITE(numout,*) ' q_i_b : ', q_i_b(jiindex_1d,nlay_i+1) 606 ! WRITE(numout,*) ' fc_bo_i : ', fc_bo_i(jiindex_1d) 607 ! WRITE(numout,*) ' fbif_1d : ', fbif_1d(jiindex_1d) 608 ! WRITE(numout,*) ' qlbbq_1d : ', qlbbq_1d(jiindex_1d) 609 ! ENDIF 610 611 !----------------------------------------- 612 ! Basal growth - salinity varying in time 613 !----------------------------------------- 451 !------------------------------------------------- 452 ! 4.1 Basal growth - (b) salinity varying in time 453 !------------------------------------------------- 614 454 IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN 615 455 ! the growth rate (dh_i_bott) is function of the new ice … … 617 457 ! salinity (snewice). snewice depends on dh_i_bott 618 458 ! it converges quickly, so, no problem 459 ! See Vancoppenolle et al., OM08 for more info on this 619 460 620 461 ! Initial value (tested 1D, can be anything between 1 and 20) … … 636 477 ( t_bo_b(ji) - rtt ) ) & 637 478 - rcp * ( ztmelts-rtt ) ) 638 ! !+++++++639 ! IF (ji.EQ.jiindex_1d) THEN640 ! WRITE(numout,*) ' ztmelts : ', ztmelts641 ! WRITE(numout,*) ' t_bo_b : ', t_bo_b(ji)642 ! WRITE(numout,*) ' q_i_b(ji,nlay_i+1) : ', q_i_b(ji,nlay_i+1)643 ! ENDIF644 ! !+++++++645 479 ! Bottom growth rate = - F*dt / q 646 480 dh_i_bott(ji) = - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) & … … 650 484 ! zswi12 (1) if dh_i_bott/rdt .LT. 3.6e-7 and .GT. 2.0e-8 651 485 ! zswi1 (1) if dh_i_bott/rdt .LT. 2.0e-8 652 ! zgrr = MAX( dh_i_bott(ji) / rdt_ice, 0.0 )653 486 zgrr = MIN( 1.0e-3, MAX ( dh_i_bott(ji) / rdt_ice , zeps ) ) 654 487 zswi2 = MAX( zzero , SIGN( zone , zgrr - 3.6e-7 ) ) 655 488 zswi12 = MAX( zzero , SIGN( zone , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 ) 656 489 zswi1 = 1. - zswi2 * zswi12 657 ! IF (ji.EQ.jiindex_1d) THEN658 ! WRITE(numout,*) ' zgrr : ', zgrr659 ! WRITE(numout,*) ' zswi2 : ', zswi2, ' zswi12 :', zswi12, ' zswi1 :', zswi1660 ! ENDIF661 490 zfracs = zswi1 * 0.12 + & 662 491 zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) ) + & … … 673 502 IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0 ) THEN 674 503 ! New ice salinity must not exceed 15 psu 675 zoldsinew = s_i_new(ji)676 504 s_i_new(ji) = MIN( s_i_new(ji), s_i_max ) 677 505 ! Metling point in K … … 684 512 - rcp * ( ztmelts-rtt ) ) 685 513 ! Basal growth rate = - F*dt / q 686 ! sometimes, growth is very high because of very high bottom cond687 ! flux... this is dangerous688 514 dh_i_bott(ji) = - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + & 689 515 qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 690 ! new lines 691 !----------------- 692 ! Salinity update 693 !----------------- 516 ! Salinity update 694 517 ! entrapment during bottom growth 695 518 dsm_i_se_1d(ji) = ( s_i_new(ji)*dh_i_bott(ji) + & … … 701 524 ENDIF ! num_sal 702 525 703 ! IF ( jiindex_1d .GT. 0 ) THEN 704 ! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(jiindex_1d) 705 ! WRITE(numout,*) ' s_i_new : ', s_i_new(jiindex_1d) 706 ! ENDIF 707 708 !------------ 709 ! Basal melt 710 !------------ 711 !++++++ 526 !---------------- 527 ! 4.2 Basal melt 528 !---------------- 712 529 meance_dh = 0.0 713 530 numce_dh = 0 714 531 innermelt(:) = 0 715 !++++++716 532 717 533 DO ji = kideb, kiut … … 729 545 END DO 730 546 731 ! IF ( jiindex_1d .GT. 0 ) THEN732 ! WRITE(numout,*) ' zqfont_bo : ', zqfont_bo(jiindex_1d) / rdt_ice733 ! WRITE(numout,*) ' fc_bo_i : ', fc_bo_i(jiindex_1d)734 ! WRITE(numout,*) ' fbif : ', fbif_1d(jiindex_1d)735 ! WRITE(numout,*) ' qlbbq : ', qlbbq_1d(jiindex_1d)736 ! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(jiindex_1d)737 ! ENDIF738 739 547 DO jk = nlay_i, 1, -1 740 548 DO ji = kideb, kiut 741 549 IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0 ) THEN 742 550 ztmelts = - tmut * s_i_b(ji,jk) + rtt 743 ! sometimes internal temperature gt melting point744 ! the whole layer is removed745 ! IF (ji.EQ.jiindex_1d) THEN746 ! WRITE(numout,*) ' jk : ', jk747 ! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)748 ! WRITE(numout,*) ' zqfont_bo : ', zqfont_bo(ji)749 ! WRITE(numout,*) ' t_i_b : ', t_i_b(ji,jk)750 ! WRITE(numout,*) ' q_i_b : ', q_i_b(ji,jk)751 ! ENDIF752 551 IF ( t_i_b(ji,jk) .GE. ztmelts ) THEN 753 552 zdeltah(ji,jk) = - zh_i(ji) 754 ! zqfont_bo(ji) = zqfont_bo(ji)755 553 dh_i_bott(ji) = dh_i_bott(ji) + zdeltah(ji,jk) 756 !++++++757 554 innermelt(ji) = 1 758 !++++++759 ! IF (ji.EQ.jiindex_1d) THEN760 ! WRITE(numout,*) ' Inner melt: ', innermelt(ji)761 ! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)762 ! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)763 ! ENDIF764 555 ELSE ! normal ablation 765 556 zdeltah(ji,jk) = - zqfont_bo(ji) / q_i_b(ji,jk) … … 774 565 zjj = ( npb(ji) - 1 ) / jpi + 1 775 566 zfsalt_melt(ji) = zfsalt_melt(ji) + & 776 ! ( sss_io(zji,zjj) - s_i_b(ji,jk) ) * &777 567 ( sss_io(zji,zjj) - sm_i_b(ji) ) * & 778 568 a_i_b(ji) * & 779 569 MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice 780 ! IF (ji.EQ.jiindex_1d) THEN781 ! WRITE(numout,*) ' No inner melt '782 ! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)783 ! WRITE(numout,*) ' zqfont_bo : ', zqfont_bo(ji)784 ! WRITE(numout,*) ' t_i_b : ', t_i_b(ji,jk)785 ! WRITE(numout,*) ' q_i_b : ', q_i_b(ji,jk)786 ! WRITE(numout,*)787 ! ENDIF788 570 ENDIF 789 ! IF ( ji .EQ. jiindex_1d ) THEN790 ! WRITE(numout,*) ' basal melt '791 ! WRITE(numout,*) ' sss_io : ', sss_io(zji,zjj)792 ! WRITE(numout,*) ' zfsalt_melt : ', zfsalt_melt(ji)793 ! WRITE(numout,*) ' zfsalt_melt good units : ', zfsalt_melt(ji) / ( sss_io(zji,zjj) + epsi16 )794 ! ENDIF795 !+++++796 571 ENDIF 797 572 END DO ! ji 798 573 END DO ! jk 799 574 575 !------------------- 576 ! Conservation test 577 !------------------- 800 578 IF ( con_i ) THEN 801 579 DO ji = kideb, kiut 802 580 IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0 ) THEN 803 !-------------------804 ! Conservation test805 !-------------------806 581 IF ( ( zfbase(ji) + zdq_i(ji) ) .GE. 1.0e-3 ) THEN 807 582 numce_dh = numce_dh + 1 808 583 meance_dh = meance_dh + zfbase(ji) + zdq_i(ji) 809 584 ENDIF 810 IF ( zfbase(ji) + zdq_i(ji) .GE. 1.0e-3 ) THEN!.GE. 1.0e-3 ) .AND. & 811 ! ( ( ht_i_b(ji) + dh_i_bott(ji) ) .GE. 1.0e-6 ) ) THEN 585 IF ( zfbase(ji) + zdq_i(ji) .GE. 1.0e-3 ) THEN 812 586 WRITE(numout,*) ' ALERTE heat loss for basal melt ' 813 587 WRITE(numout,*) ' zji, zjj, jl :', zji, zjj, jl … … 834 608 ENDIF ! con_i 835 609 836 ! IF ( jiindex_1d .GT. 0 ) THEN837 ! WRITE(numout,*) ' zqfont_bo : ', zqfont_bo(jiindex_1d) / rdt_ice838 ! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(jiindex_1d)839 ! ENDIF840 610 ! 841 611 !------------------------------------------------------------------------------! … … 843 613 !------------------------------------------------------------------------------! 844 614 ! 845 !------------------------------------------ 846 ! Excessive ablation in a 1-category model847 !------------------------------------------ 615 !---------------------------------------------- 616 ! 5.1 Excessive ablation in a 1-category model 617 !---------------------------------------------- 848 618 849 619 DO ji = kideb, kiut … … 872 642 END DO 873 643 874 ! IF (jiindex_1d .GT. 0 ) THEN 875 ! WRITE(numout,*) ' ht_i_b : ', ht_i_b(jiindex_1d) 876 ! WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d) 877 ! ENDIF 878 879 !-------------------------- 880 ! If too much ice melts 881 !-------------------------- 644 !----------------------------------- 645 ! 5.2 More than available ice melts 646 !----------------------------------- 882 647 ! then heat applied minus heat content at previous time step 883 648 ! should equal heat remaining … … 890 655 zhgnew(ji) = zihgnew * zhgnew(ji) ! ice thickness is put to 0 891 656 ! remaining heat 892 zfdt_final(ji) = ( 1.0 - zihgnew ) * ( zqfont_su(ji) + zqfont_bo(ji) ) ! & 893 ! + zihgnew * MAX ( zfdt_init(ji) - zqt_i(ji) - zqt_s(ji), 0.0 ) 894 895 ! !++++++++++ 896 ! IF (ji.EQ.jiindex_1d) THEN 897 ! WRITE(numout,*) ' zfdt_init : ', zfdt_init(ji) / rdt_ice 898 ! WRITE(numout,*) ' zqt_i : ', zqt_i(ji) / rdt_ice 899 ! WRITE(numout,*) ' 1. zqt_s : ', zqt_s(ji) / rdt_ice 900 ! WRITE(numout,*) ' zqt : ', ( zqt_i(ji) + zqt_s(ji) ) / rdt_ice 901 ! WRITE(numout,*) ' zhgnew : ', zhgnew(ji) 902 ! WRITE(numout,*) ' zihgnew : ', zihgnew 903 ! WRITE(numout,*) ' dh_i_surf : ', dh_i_surf(ji) 904 ! WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji) 905 ! WRITE(numout,*) ' ht_s_b : ', ht_s_b(ji) 906 ! WRITE(numout,*) ' 1. zfdt_final: ', zfdt_final(ji) / rdt_ice 907 ! WRITE(numout,*) ' zqfont_su : ', zqfont_su(jiindex_1d) / rdt_ice 908 ! WRITE(numout,*) ' zqfont_bo : ', zqfont_bo(jiindex_1d) / rdt_ice 909 ! ENDIF 910 ! !++++++++++ 657 zfdt_final(ji) = ( 1.0 - zihgnew ) * ( zqfont_su(ji) + zqfont_bo(ji) ) 911 658 912 659 ! If snow remains, energy is used to melt snow … … 923 670 zqt_s(ji) = zqt_s(ji) * ht_s_b(ji) 924 671 925 ! IF (ji.EQ.jiindex_1d) THEN926 ! WRITE(numout,*) ' 2. zfdt_final : ', zfdt_final(ji) / rdt_ice927 ! WRITE(numout,*) ' ht_s_b : ', ht_s_b(ji)928 ! WRITE(numout,*) ' zhgnew : ', zhgnew(ji)929 ! WRITE(numout,*) ' 2. zqt_s : ', zqt_s(ji) / rdt_ice930 ! WRITE(numout,*) ' zdhnm : ', zdhnm931 ! WRITE(numout,*)932 ! ENDIF933 934 672 ! Mass variations of ice and snow 935 673 !--------------------------------- … … 942 680 ! Remaining heat to the ocean 943 681 !--------------------------------- 944 ! check the switches here945 682 ! focea is in W.m-2 * dt 946 947 683 focea(ji) = - zfdt_final(ji) / rdt_ice 948 684 949 ! IF ( (zji.eq.jiindex) .AND. (zjj.eq.jjindex) ) THEN 950 ! WRITE(numout,*) ' focea : ', focea(ji) 951 ! ENDIF 952 953 END DO 954 955 ! IF (jiindex_1d .GT. 0 ) THEN 956 ! WRITE(numout,*) ' ht_i_b : ', ht_i_b(jiindex_1d) 957 ! WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d) 958 ! ENDIF 685 END DO 959 686 960 687 ftotal_fin (:) = zfdt_final(:) / rdt_ice 961 ! IF (jiindex_1d .GT. 1 ) WRITE(numout,*) ' ftotal_fin : ', ftotal_fin(jiindex_1d)962 688 963 689 !--------------------------- … … 986 712 focea(ji) * a_i_b(ji) * rdt_ice 987 713 988 ! ! astarojna989 714 zihic = 1.0 - MAX( zzero , SIGN( zone , -ht_i_b(ji) ) ) 990 715 ! equals 0 if ht_i = 0, 1 if ht_i gt 0 991 ! fstbif_1d est cumulatif merde992 716 fscbq_1d(ji) = a_i_b(ji) * fstbif_1d(ji) 993 ! here there is a bug994 717 qldif_1d(ji) = qldif_1d(ji) & 995 718 + fsup(ji) + ( 1.0 - zihgnew ) * focea(ji) * a_i_b(ji) & … … 1011 734 t_i_b(ji,jk) = zihgnew * t_i_b(ji,jk) + ( 1.0 - zihgnew ) * rtt 1012 735 q_i_b(ji,jk) = zihgnew * q_i_b(ji,jk) 1013 ! IF (ji.eq.jiindex_1d) THEN1014 ! WRITE(numout,*) ' jk : ', jk1015 ! WRITE(numout,*) ' q_i_b : ', q_i_b(ji,jk)1016 ! WRITE(numout,*) ' t_i_b : ', t_i_b(jiindex_1d, 1:nlay_i)1017 ! WRITE(numout,*) ' zihgnew : ', zihgnew, ' zhgnew : ', &1018 ! zhgnew(ji)1019 ! ENDIF1020 736 END DO 1021 737 END DO ! ji 1022 738 1023 ! IF (jiindex_1d.GT.0) THEN1024 ! WRITE(numout,*) ' --- Vertical profile : '1025 ! WRITE(numout,*) ' q_i_b : ', q_i_b(jiindex_1d, 1:nlay_i)1026 ! ENDIF1027 1028 739 DO ji = kideb, kiut 1029 740 ht_i_b(ji) = zhgnew(ji) 1030 741 END DO ! ji 1031 1032 ! IF (jiindex_1d .GT. 0 ) THEN1033 ! WRITE(numout,*) ' ht_i_b : ', ht_i_b(jiindex_1d)1034 ! WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d)1035 ! ENDIF1036 742 ! 1037 743 !------------------------------------------------------------------------------| … … 1065 771 zji = MOD( npb(ji) - 1, jpi ) + 1 1066 772 zjj = ( npb(ji) - 1 ) / jpi + 1 1067 ! zsm_snowice = MIN ( s_i_max, ( rhoic - rhosn ) / rhoic * &1068 ! sss_io(zji,zjj) )1069 773 1070 774 zsm_snowice = ( rhoic - rhosn ) / rhoic * & … … 1115 819 END DO !ji 1116 820 1117 ! IF (jiindex_1d .GT. 1 ) THEN1118 ! WRITE(numout,*) ' ht_i_b : ', ht_i_b(jiindex_1d)1119 ! WRITE(numout,*) ' ht_s_b : ', ht_s_b(jiindex_1d)1120 ! ENDIF1121 1122 821 END SUBROUTINE lim_thd_dh 1123 822 #else
Note: See TracChangeset
for help on using the changeset viewer.
