- Timestamp:
- 2013-11-07T11:01:27+01:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r3808 r4161 6 6 !! History : LIM ! 2003-05 (M. Vancoppenolle) Original code in 1D 7 7 !! ! 2005-06 (M. Vancoppenolle) 3D version 8 !! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in rdm snif & rdmicif8 !! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in rdm_snw & rdm_ice 9 9 !! 3.4 ! 2011-02 (G. Madec) dynamical allocation 10 10 !! 3.5 ! 2012-10 (G. Madec & co) salt flux + bug fixes … … 39 39 40 40 !!---------------------------------------------------------------------- 41 !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011)41 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 42 42 !! $Id$ 43 43 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 73 73 !! 74 74 INTEGER :: ji , jk ! dummy loop indices 75 INTEGER :: ii, ij ! 2D corresponding indices to ji75 INTEGER :: ii, ij ! 2D corresponding indices to ji 76 76 INTEGER :: isnow ! switch for presence (1) or absence (0) of snow 77 77 INTEGER :: isnowic ! snow ice formation not … … 84 84 REAL(wp) :: zdhnm, zhnnew, zhisn, zihic, zzc ! 85 85 REAL(wp) :: zfracs ! fractionation coefficient for bottom salt entrapment 86 REAL(wp) :: zds ! increment of bottom ice salinity87 86 REAL(wp) :: zcoeff ! dummy argument for snowfall partitioning over ice and leads 88 87 REAL(wp) :: zsm_snowice ! snow-ice salinity … … 107 106 REAL(wp), POINTER, DIMENSION(:) :: zdh_s_pre ! snow precipitation 108 107 REAL(wp), POINTER, DIMENSION(:) :: zdh_s_sub ! snow sublimation 109 REAL(wp), POINTER, DIMENSION(:) :: zsfx_melt ! salt flux due to ice melt110 108 111 109 REAL(wp), POINTER, DIMENSION(:,:) :: zdeltah … … 120 118 REAL(wp), POINTER, DIMENSION(:,:) :: zqt_i_lay ! total ice heat content 121 119 120 ! mass and salt flux (clem) 121 REAL(wp) :: zdvres, zdvsur, zdvbot 122 REAL(wp), POINTER, DIMENSION(:) :: zviold, zvsold ! old ice volume... 123 122 124 ! Heat conservation 123 125 INTEGER :: num_iter_max, numce_dh … … 128 130 129 131 CALL wrk_alloc( jpij, zh_i, zh_s, ztfs, zhsold, zqprec, zqfont_su, zqfont_bo, z_f_surf, zhgnew, zfmass_i ) 130 CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, z sfx_melt, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy )132 CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy ) 131 133 CALL wrk_alloc( jpij, zinnermelt, zfbase, zdq_i ) 132 134 CALL wrk_alloc( jpij, jkmax, zdeltah, zqt_i_lay ) 133 135 134 zsfx_melt (:) = 0._wp 136 CALL wrk_alloc( jpij, zviold, zvsold ) ! clem 137 135 138 ftotal_fin(:) = 0._wp 136 139 zfdt_init (:) = 0._wp 137 140 zfdt_final(:) = 0._wp 138 141 142 dh_i_surf (:) = 0._wp 143 dh_i_bott (:) = 0._wp 144 dh_snowice(:) = 0._wp 145 139 146 DO ji = kideb, kiut 140 147 old_ht_i_b(ji) = ht_i_b(ji) 141 148 old_ht_s_b(ji) = ht_s_b(ji) 149 zviold(ji) = a_i_b(ji) * ht_i_b(ji) ! clem 150 zvsold(ji) = a_i_b(ji) * ht_s_b(ji) ! clem 142 151 END DO 143 152 ! … … 164 173 ! 165 174 DO ji = kideb, kiut ! Layer thickness 166 zh_i(ji) = ht_i_b(ji) / nlay_i167 zh_s(ji) = ht_s_b(ji) / nlay_s175 zh_i(ji) = ht_i_b(ji) / REAL( nlay_i ) 176 zh_s(ji) = ht_s_b(ji) / REAL( nlay_s ) 168 177 END DO 169 178 ! … … 171 180 DO jk = 1, nlay_s 172 181 DO ji = kideb, kiut 173 zqt_s(ji) = zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s182 zqt_s(ji) = zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / REAL( nlay_s ) 174 183 END DO 175 184 END DO … … 178 187 DO jk = 1, nlay_i 179 188 DO ji = kideb, kiut 180 zzc = q_i_b(ji,jk) * ht_i_b(ji) / nlay_i189 zzc = q_i_b(ji,jk) * ht_i_b(ji) / REAL( nlay_i ) 181 190 zqt_i(ji) = zqt_i(ji) + zzc 182 191 zqt_i_lay(ji,jk) = zzc … … 244 253 zhn = 1.0 - MAX( zzero , SIGN( zone , - zhsnew ) ) 245 254 ht_s_b(ji) = MAX( zzero , zhsnew ) 255 ! we recompute dh_s_tot (clem) 256 dh_s_tot (ji) = ht_s_b(ji) - zhsold(ji) 246 257 ! Volume and mass variations of snow 247 258 dvsbq_1d (ji) = a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji) - zdh_s_pre(ji) ) 248 259 dvsbq_1d (ji) = MIN( zzero, dvsbq_1d(ji) ) 249 rdm_snw_1d(ji) = rdm_snw_1d(ji) + rhosn * dvsbq_1d(ji)260 !clem rdm_snw_1d(ji) = rdm_snw_1d(ji) + rhosn * dvsbq_1d(ji) 250 261 END DO ! ji 251 262 … … 254 265 !-------------------------- 255 266 DO ji = kideb, kiut 256 dh_i_surf(ji) = 0._wp257 267 z_f_surf (ji) = zqfont_su(ji) * r1_rdtice ! heat conservation test 258 268 zdq_i (ji) = 0._wp … … 272 282 zdq_i (ji ) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) * r1_rdtice 273 283 ! 274 ! ! contribution to ice-ocean salt flux 275 zsfx_melt(ji) = zsfx_melt(ji) - sm_i_b(ji) * a_i_b(ji) * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic * r1_rdtice 284 ! clem 285 sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) & 286 & * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic / rdt_ice 276 287 END DO 277 288 END DO … … 334 345 DO ji = kideb,kiut 335 346 q_s_b (ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) 336 zqt_dummy(ji) = zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s! heat conservation347 zqt_dummy(ji) = zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / REAL( nlay_s ) ! heat conservation 337 348 END DO 338 349 END DO … … 375 386 ! Basal growth rate = - F*dt / q 376 387 dh_i_bott(ji) = - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 388 sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * dh_i_bott(ji) * rhoic * r1_rdtice 377 389 ENDIF 378 390 END DO … … 416 428 zfracs = zswi1 * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) ) & 417 429 & + zswi2 * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) 418 z ds = zfracs * sss_m(ii,ij) - s_i_new(ji)430 zfracs = MIN( 0.5 , zfracs ) 419 431 s_i_new(ji) = zfracs * sss_m(ii,ij) 420 432 ENDIF ! fc_bo_i … … 425 437 DO ji = kideb, kiut 426 438 IF( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0 ) THEN 427 ! New ice salinity must not exceed 15psu439 ! New ice salinity must not exceed 20 psu 428 440 s_i_new(ji) = MIN( s_i_new(ji), s_i_max ) 429 441 ! Metling point in K … … 437 449 ! Salinity update 438 450 ! entrapment during bottom growth 439 dsm_i_se_1d(ji) = ( s_i_new(ji) * dh_i_bott(ji) + sm_i_b(ji) * ht_i_b(ji) ) & 440 & / MAX( ht_i_b(ji) + dh_i_bott(ji) ,epsi13 ) - sm_i_b(ji) 451 sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * dh_i_bott(ji) * rhoic * r1_rdtice 441 452 ENDIF ! heat budget 442 453 END DO … … 476 487 zdq_i (ji ) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) * r1_rdtice 477 488 ENDIF 478 ! contribution to salt flux 479 zsfx_melt(ji) = zsfx_melt(ji) - sm_i_b(ji) * a_i_b(ji) * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic * r1_rdtice 489 ! clem: contribution to salt flux 490 sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) & 491 & * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic * r1_rdtice 480 492 ENDIF 481 493 END DO ! ji … … 528 540 ELSE ; zdhbf = dh_i_bott(ji) 529 541 ENDIF 542 zdvres = zdhbf - dh_i_bott(ji) 543 dh_i_bott(ji) = zdhbf 544 sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvres * rhoic * r1_rdtice 530 545 ! ! excessive energy is sent to lateral ablation 531 fsup (ji) = rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi13 ) & 532 & * ( zdhbf - dh_i_bott(ji) ) * r1_rdtice 533 dh_i_bott(ji) = zdhbf 534 ! !since ice volume is only used for outputs, we keep it global for all categories 535 dvbbq_1d (ji) = a_i_b(ji) * dh_i_bott(ji) 536 ! !new ice thickness 537 zhgnew (ji) = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji) 538 ! ! diagnostic ( bottom ice growth ) 539 ii = MOD( npb(ji) - 1, jpi ) + 1 540 ij = ( npb(ji) - 1 ) / jpi + 1 541 diag_bot_gr(ii,ij) = diag_bot_gr(ii,ij) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice 542 diag_sur_me(ii,ij) = diag_sur_me(ii,ij) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) * r1_rdtice 543 diag_bot_me(ii,ij) = diag_bot_me(ii,ij) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice 546 fsup (ji) = rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi13 ) * zdvres * r1_rdtice 544 547 END DO 545 548 … … 552 555 ! Adapt the remaining energy if too much ice melts 553 556 !-------------------------------------------------- 554 zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) ! =1 if ice 555 ! 0 if no more ice 556 zhgnew (ji) = zihgnew * zhgnew(ji) ! ice thickness is put to 0 557 ! remaining heat 557 zdvres = MAX( 0._wp, - ht_i_b(ji) - dh_i_surf(ji) - dh_i_bott(ji) ) 558 zdvsur = MIN( 0._wp, dh_i_surf(ji) + zdvres ) - dh_i_surf(ji) ! fill the surface first 559 zdvbot = MAX( 0._wp, zdvres - zdvsur ) ! then the bottom 560 dh_i_surf (ji) = dh_i_surf(ji) + zdvsur ! clem 561 dh_i_bott (ji) = dh_i_bott(ji) + zdvbot ! clem 562 563 ! new ice thickness (clem) 564 zhgnew(ji) = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji) 565 zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice 566 zhgnew(ji) = zihgnew * zhgnew(ji) ! ice thickness is put to 0 567 568 ! !since ice volume is only used for outputs, we keep it global for all categories 569 dvbbq_1d (ji) = a_i_b(ji) * dh_i_bott(ji) 570 571 ! remaining heat 558 572 zfdt_final(ji) = ( 1.0 - zihgnew ) * ( zqfont_su(ji) + zqfont_bo(ji) ) 559 573 … … 569 583 ht_s_b(ji) = MAX( zzero , zhnfi ) 570 584 zqt_s(ji) = zqt_s(ji) * ht_s_b(ji) 585 ! we recompute dh_s_tot (clem) 586 dh_s_tot (ji) = ht_s_b(ji) - zhsold(ji) 571 587 572 588 ! Mass variations of ice and snow … … 579 595 ! 580 596 ! ! mass variation cumulated over category 581 rdm_snw_1d(ji) = rdm_snw_1d(ji) + zzfmass_s ! snow582 rdm_ice_1d(ji) = rdm_ice_1d(ji) + zzfmass_i ! ice597 !clem rdm_snw_1d(ji) = rdm_snw_1d(ji) + zzfmass_s ! snow 598 !clem rdm_ice_1d(ji) = rdm_ice_1d(ji) + zzfmass_i ! ice 583 599 584 600 ! Remaining heat to the ocean … … 586 602 focea(ji) = - zfdt_final(ji) * r1_rdtice ! focea is in W.m-2 * dt 587 603 604 ! residual salt flux (clem) 605 !-------------------------- 606 ! surface 607 sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvsur * rhoic * r1_rdtice 608 ! bottom 609 IF ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) >= 0._wp ) THEN ! melting 610 sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvbot * rhoic * r1_rdtice 611 ELSE ! growth 612 sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * zdvbot * rhoic * r1_rdtice 613 ENDIF 614 ! 615 ! diagnostic ( bottom ice growth ) 616 ii = MOD( npb(ji) - 1, jpi ) + 1 617 ij = ( npb(ji) - 1 ) / jpi + 1 618 diag_bot_gr(ii,ij) = diag_bot_gr(ii,ij) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice 619 diag_sur_me(ii,ij) = diag_sur_me(ii,ij) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) * r1_rdtice 620 diag_bot_me(ii,ij) = diag_bot_me(ii,ij) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice 588 621 END DO 589 622 … … 591 624 592 625 !--------------------------- 593 ! Salt flux andheat fluxes626 ! heat fluxes 594 627 !--------------------------- 595 628 DO ji = kideb, kiut 596 629 zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) ! =1 if ice 597 !598 ! Salt flux599 sfx_thd_1d(ji) = sfx_thd_1d(ji) + zihgnew * zsfx_melt(ji) &600 & - (1.0 - zihgnew) * zfmass_i (ji) * sm_i_b(ji) * r1_rdtice601 630 ! 602 631 ! Heat flux … … 646 675 dmgwi_1d (ji) = dmgwi_1d(ji) + a_i_b(ji) * ( ht_s_b(ji) - zhnnew ) * rhosn 647 676 648 ! All snow is thrown in the ocean, and seawater is taken to replace the volume 649 rdm_ice_1d(ji) = rdm_ice_1d(ji) + a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic * ( 1. - rhosn / rhoic ) 650 rdm_snw_1d(ji) = rdm_snw_1d(ji) + a_i_b(ji) * ( zhnnew - ht_s_b(ji) ) * rhosn 677 !clem rdm_ice_1d(ji) = rdm_ice_1d(ji) + a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic 678 !clem rdm_snw_1d(ji) = rdm_snw_1d(ji) + a_i_b(ji) * ( zhnnew - ht_s_b(ji) ) * rhosn 651 679 652 680 ! Equivalent salt flux (1) Snow-ice formation component … … 658 686 ELSE ; zsm_snowice = sm_i_b(ji) 659 687 ENDIF 660 sfx_thd_1d(ji) = sfx_thd_1d(ji) - zsm_snowice * a_i_b(ji) * dh_snowice(ji) * rhoic * r1_rdtice661 !662 688 ! entrapment during snow ice formation 663 i_ice_switch = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - ht_i_b(ji) + 1.0e-6 ) ) 664 isnowic = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - dh_snowice(ji) ) ) * i_ice_switch 665 IF( num_sal == 2 ) & 666 dsm_i_si_1d(ji) = ( zsm_snowice * dh_snowice(ji) & 667 & + sm_i_b(ji) * ht_i_b(ji) / MAX( ht_i_b(ji) + dh_snowice(ji), epsi13 ) & 668 & - sm_i_b(ji) ) * isnowic 689 ! clem: new salinity difference stored (to be used in limthd_ent.F90) 690 IF ( num_sal == 2 ) THEN 691 i_ice_switch = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - zhgnew(ji) + epsi13 ) ) 692 ! salinity dif due to snow-ice formation 693 dsm_i_si_1d(ji) = ( zsm_snowice - sm_i_b(ji) ) * dh_snowice(ji) / MAX( zhgnew(ji), epsi13 ) * i_ice_switch 694 ! salinity dif due to bottom growth 695 IF ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) < 0._wp ) THEN 696 dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_b(ji) ) * dh_i_bott(ji) / MAX( zhgnew(ji), epsi13 ) * i_ice_switch 697 ENDIF 698 ENDIF 669 699 670 700 ! Actualize new snow and ice thickness. … … 680 710 diag_sni_gr(ii,ij) = diag_sni_gr(ii,ij) + dh_snowice(ji)*a_i_b(ji) * r1_rdtice 681 711 ! 712 ! salt flux 713 sfx_thd_1d(ji) = sfx_thd_1d(ji) - zsm_snowice * a_i_b(ji) * dh_snowice(ji) * rhoic * r1_rdtice 714 !-------------------------------- 715 ! Update mass fluxes (clem) 716 !-------------------------------- 717 rdm_ice_1d(ji) = rdm_ice_1d(ji) + ( a_i_b(ji) * ht_i_b(ji) - zviold(ji) ) * rhoic 718 rdm_snw_1d(ji) = rdm_snw_1d(ji) + ( a_i_b(ji) * ht_s_b(ji) - zvsold(ji) ) * rhosn 719 682 720 END DO !ji 683 721 ! 684 722 CALL wrk_dealloc( jpij, zh_i, zh_s, ztfs, zhsold, zqprec, zqfont_su, zqfont_bo, z_f_surf, zhgnew, zfmass_i ) 685 CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, z sfx_melt, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy )723 CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy ) 686 724 CALL wrk_dealloc( jpij, zinnermelt, zfbase, zdq_i ) 687 725 CALL wrk_dealloc( jpij, jkmax, zdeltah, zqt_i_lay ) 726 ! 727 CALL wrk_dealloc( jpij, zviold, zvsold ) ! clem 688 728 ! 689 729 END SUBROUTINE lim_thd_dh
Note: See TracChangeset
for help on using the changeset viewer.
