- Timestamp:
- 2013-09-25T16:38:24+02:00 (11 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r3808 r4045 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 … … 107 107 REAL(wp), POINTER, DIMENSION(:) :: zdh_s_pre ! snow precipitation 108 108 REAL(wp), POINTER, DIMENSION(:) :: zdh_s_sub ! snow sublimation 109 REAL(wp), POINTER, DIMENSION(:) :: zsfx_melt ! salt flux due to ice melt110 109 111 110 REAL(wp), POINTER, DIMENSION(:,:) :: zdeltah … … 120 119 REAL(wp), POINTER, DIMENSION(:,:) :: zqt_i_lay ! total ice heat content 121 120 121 ! mass and salt flux (clem) 122 REAL(wp) :: zdvres, zdvsur, zdvbot 123 REAL(wp), POINTER, DIMENSION(:) :: zviold, zvsold ! old ice volume... 124 122 125 ! Heat conservation 123 126 INTEGER :: num_iter_max, numce_dh … … 128 131 129 132 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 )133 CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy ) 131 134 CALL wrk_alloc( jpij, zinnermelt, zfbase, zdq_i ) 132 135 CALL wrk_alloc( jpij, jkmax, zdeltah, zqt_i_lay ) 133 136 134 zsfx_melt (:) = 0._wp 137 CALL wrk_alloc( jpij, zviold, zvsold ) ! clem 138 135 139 ftotal_fin(:) = 0._wp 136 140 zfdt_init (:) = 0._wp 137 141 zfdt_final(:) = 0._wp 138 142 143 dh_i_surf (:) = 0._wp 144 dh_i_bott (:) = 0._wp 145 dh_snowice(:) = 0._wp 146 139 147 DO ji = kideb, kiut 140 148 old_ht_i_b(ji) = ht_i_b(ji) 141 149 old_ht_s_b(ji) = ht_s_b(ji) 150 zviold(ji) = a_i_b(ji) * ht_i_b(ji) ! clem 151 zvsold(ji) = a_i_b(ji) * ht_s_b(ji) ! clem 142 152 END DO 143 153 ! … … 164 174 ! 165 175 DO ji = kideb, kiut ! Layer thickness 166 zh_i(ji) = ht_i_b(ji) / nlay_i167 zh_s(ji) = ht_s_b(ji) / nlay_s176 zh_i(ji) = ht_i_b(ji) / REAL( nlay_i ) 177 zh_s(ji) = ht_s_b(ji) / REAL( nlay_s ) 168 178 END DO 169 179 ! … … 171 181 DO jk = 1, nlay_s 172 182 DO ji = kideb, kiut 173 zqt_s(ji) = zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s183 zqt_s(ji) = zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / REAL( nlay_s ) 174 184 END DO 175 185 END DO … … 178 188 DO jk = 1, nlay_i 179 189 DO ji = kideb, kiut 180 zzc = q_i_b(ji,jk) * ht_i_b(ji) / nlay_i190 zzc = q_i_b(ji,jk) * ht_i_b(ji) / REAL( nlay_i ) 181 191 zqt_i(ji) = zqt_i(ji) + zzc 182 192 zqt_i_lay(ji,jk) = zzc … … 244 254 zhn = 1.0 - MAX( zzero , SIGN( zone , - zhsnew ) ) 245 255 ht_s_b(ji) = MAX( zzero , zhsnew ) 256 ! we recompute dh_s_tot (clem) 257 dh_s_tot (ji) = ht_s_b(ji) - zhsold(ji) 246 258 ! Volume and mass variations of snow 247 259 dvsbq_1d (ji) = a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji) - zdh_s_pre(ji) ) 248 260 dvsbq_1d (ji) = MIN( zzero, dvsbq_1d(ji) ) 249 rdm_snw_1d(ji) = rdm_snw_1d(ji) + rhosn * dvsbq_1d(ji)261 !clem rdm_snw_1d(ji) = rdm_snw_1d(ji) + rhosn * dvsbq_1d(ji) 250 262 END DO ! ji 251 263 … … 254 266 !-------------------------- 255 267 DO ji = kideb, kiut 256 dh_i_surf(ji) = 0._wp257 268 z_f_surf (ji) = zqfont_su(ji) * r1_rdtice ! heat conservation test 258 269 zdq_i (ji) = 0._wp … … 272 283 zdq_i (ji ) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) * r1_rdtice 273 284 ! 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 285 ! clem 286 sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) & 287 & * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic / rdt_ice 276 288 END DO 277 289 END DO … … 334 346 DO ji = kideb,kiut 335 347 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 conservation348 zqt_dummy(ji) = zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / REAL( nlay_s ) ! heat conservation 337 349 END DO 338 350 END DO … … 375 387 ! Basal growth rate = - F*dt / q 376 388 dh_i_bott(ji) = - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 389 sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * dh_i_bott(ji) * rhoic * r1_rdtice 377 390 ENDIF 378 391 END DO … … 437 450 ! Salinity update 438 451 ! 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) 452 sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * dh_i_bott(ji) * rhoic * r1_rdtice 441 453 ENDIF ! heat budget 442 454 END DO … … 476 488 zdq_i (ji ) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) * r1_rdtice 477 489 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 490 ! clem: contribution to salt flux 491 sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) & 492 & * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic * r1_rdtice 480 493 ENDIF 481 494 END DO ! ji … … 528 541 ELSE ; zdhbf = dh_i_bott(ji) 529 542 ENDIF 543 zdvres = zdhbf - dh_i_bott(ji) 544 dh_i_bott(ji) = zdhbf 545 sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvres * rhoic * r1_rdtice 530 546 ! ! 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 547 fsup (ji) = rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi13 ) * zdvres * r1_rdtice 544 548 END DO 545 549 … … 552 556 ! Adapt the remaining energy if too much ice melts 553 557 !-------------------------------------------------- 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 558 zdvres = MAX( 0._wp, - ht_i_b(ji) - dh_i_surf(ji) - dh_i_bott(ji) ) 559 zdvsur = MIN( 0._wp, dh_i_surf(ji) + zdvres ) - dh_i_surf(ji) ! fill the surface first 560 zdvbot = MAX( 0._wp, zdvres - zdvsur ) ! then the bottom 561 dh_i_surf (ji) = dh_i_surf(ji) + zdvsur ! clem 562 dh_i_bott (ji) = dh_i_bott(ji) + zdvbot ! clem 563 564 ! new ice thickness (clem) 565 zhgnew(ji) = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji) 566 zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice 567 zhgnew(ji) = zihgnew * zhgnew(ji) ! ice thickness is put to 0 568 569 ! !since ice volume is only used for outputs, we keep it global for all categories 570 dvbbq_1d (ji) = a_i_b(ji) * dh_i_bott(ji) 571 572 ! remaining heat 558 573 zfdt_final(ji) = ( 1.0 - zihgnew ) * ( zqfont_su(ji) + zqfont_bo(ji) ) 559 574 … … 569 584 ht_s_b(ji) = MAX( zzero , zhnfi ) 570 585 zqt_s(ji) = zqt_s(ji) * ht_s_b(ji) 586 ! we recompute dh_s_tot (clem) 587 dh_s_tot (ji) = ht_s_b(ji) - zhsold(ji) 571 588 572 589 ! Mass variations of ice and snow … … 579 596 ! 580 597 ! ! 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 ! ice598 !clem rdm_snw_1d(ji) = rdm_snw_1d(ji) + zzfmass_s ! snow 599 !clem rdm_ice_1d(ji) = rdm_ice_1d(ji) + zzfmass_i ! ice 583 600 584 601 ! Remaining heat to the ocean … … 586 603 focea(ji) = - zfdt_final(ji) * r1_rdtice ! focea is in W.m-2 * dt 587 604 605 ! residual salt flux (clem) 606 !-------------------------- 607 ! surface 608 sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvsur * rhoic * r1_rdtice 609 ! bottom 610 IF ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) >= 0._wp ) THEN ! melting 611 sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvbot * rhoic * r1_rdtice 612 ELSE ! growth 613 sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * zdvbot * rhoic * r1_rdtice 614 ENDIF 615 ! 616 ! diagnostic ( bottom ice growth ) 617 ii = MOD( npb(ji) - 1, jpi ) + 1 618 ij = ( npb(ji) - 1 ) / jpi + 1 619 diag_bot_gr(ii,ij) = diag_bot_gr(ii,ij) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice 620 diag_sur_me(ii,ij) = diag_sur_me(ii,ij) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) * r1_rdtice 621 diag_bot_me(ii,ij) = diag_bot_me(ii,ij) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice 588 622 END DO 589 623 … … 591 625 592 626 !--------------------------- 593 ! Salt flux andheat fluxes627 ! heat fluxes 594 628 !--------------------------- 595 629 DO ji = kideb, kiut 596 630 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 631 ! 602 632 ! Heat flux … … 646 676 dmgwi_1d (ji) = dmgwi_1d(ji) + a_i_b(ji) * ( ht_s_b(ji) - zhnnew ) * rhosn 647 677 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 678 !clem rdm_ice_1d(ji) = rdm_ice_1d(ji) + a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic 679 !clem rdm_snw_1d(ji) = rdm_snw_1d(ji) + a_i_b(ji) * ( zhnnew - ht_s_b(ji) ) * rhosn 651 680 652 681 ! Equivalent salt flux (1) Snow-ice formation component … … 658 687 ELSE ; zsm_snowice = sm_i_b(ji) 659 688 ENDIF 660 sfx_thd_1d(ji) = sfx_thd_1d(ji) - zsm_snowice * a_i_b(ji) * dh_snowice(ji) * rhoic * r1_rdtice661 !662 689 ! 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 690 ! clem: new salinity difference stored (to be used in limthd_ent.F90) 691 IF ( num_sal == 2 ) THEN 692 i_ice_switch = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - zhgnew(ji) + epsi13 ) ) 693 ! salinity dif due to snow-ice formation 694 dsm_i_si_1d(ji) = ( zsm_snowice - sm_i_b(ji) ) * dh_snowice(ji) / MAX( zhgnew(ji), epsi13 ) * i_ice_switch 695 ! salinity dif due to bottom growth 696 IF ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) < 0._wp ) THEN 697 dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_b(ji) ) * dh_i_bott(ji) / MAX( zhgnew(ji), epsi13 ) * i_ice_switch 698 ENDIF 699 ENDIF 669 700 670 701 ! Actualize new snow and ice thickness. … … 680 711 diag_sni_gr(ii,ij) = diag_sni_gr(ii,ij) + dh_snowice(ji)*a_i_b(ji) * r1_rdtice 681 712 ! 713 ! salt flux 714 sfx_thd_1d(ji) = sfx_thd_1d(ji) - zsm_snowice * a_i_b(ji) * dh_snowice(ji) * rhoic * r1_rdtice 715 !-------------------------------- 716 ! Update mass fluxes (clem) 717 !-------------------------------- 718 rdm_ice_1d(ji) = rdm_ice_1d(ji) + ( a_i_b(ji) * ht_i_b(ji) - zviold(ji) ) * rhoic 719 rdm_snw_1d(ji) = rdm_snw_1d(ji) + ( a_i_b(ji) * ht_s_b(ji) - zvsold(ji) ) * rhosn 720 682 721 END DO !ji 683 722 ! 684 723 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 )724 CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy ) 686 725 CALL wrk_dealloc( jpij, zinnermelt, zfbase, zdq_i ) 687 726 CALL wrk_dealloc( jpij, jkmax, zdeltah, zqt_i_lay ) 727 ! 728 CALL wrk_dealloc( jpij, zviold, zvsold ) ! clem 688 729 ! 689 730 END SUBROUTINE lim_thd_dh
Note: See TracChangeset
for help on using the changeset viewer.
