- Timestamp:
- 2015-02-05T18:54:24+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r5055 r5064 109 109 ! 1.2) Heat content 110 110 !-------------------- 111 ! Change the units of heat content; from global units to J.m3111 ! Change the units of heat content; from J/m2 to J/m3 112 112 DO jl = 1, jpl 113 113 DO jk = 1, nlay_i … … 115 115 DO ji = 1, jpi 116 116 !0 if no ice and 1 if yes 117 rswitch = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi 10 ) )117 rswitch = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi20 ) ) 118 118 !Energy of melting q(S,T) [J.m-3] 119 e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi10 ) ) * REAL( nlay_i ) 120 !convert units ! very important that this line is here 121 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac 119 e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL( nlay_i ) 122 120 END DO 123 121 END DO … … 127 125 DO ji = 1, jpi 128 126 !0 if no ice and 1 if yes 129 rswitch = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_s(ji,jj,jl) + epsi 10 ) )127 rswitch = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_s(ji,jj,jl) + epsi20 ) ) 130 128 !Energy of melting q(S,T) [J.m-3] 131 e_s(ji,jj,jk,jl) = rswitch * e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL( nlay_s ) 132 !convert units ! very important that this line is here 133 e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * unit_fac 129 e_s(ji,jj,jk,jl) = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL( nlay_s ) 134 130 END DO 135 131 END DO … … 453 449 ! Ice heat content 454 450 !------------------------ 455 ! Enthalpies are global variables we have to readjust the units (heat content in J oules)451 ! Enthalpies are global variables we have to readjust the units (heat content in J/m2) 456 452 DO jl = 1, jpl 457 453 DO jk = 1, nlay_i 458 e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * a rea(:,:) * a_i(:,:,jl) * ht_i(:,:,jl) / ( unit_fac * REAL( nlay_i ))454 e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * a_i(:,:,jl) * ht_i(:,:,jl) / REAL( nlay_i ) 459 455 END DO 460 456 END DO … … 463 459 ! Snow heat content 464 460 !------------------------ 465 ! Enthalpies are global variables we have to readjust the units (heat content in J oules)461 ! Enthalpies are global variables we have to readjust the units (heat content in J/m2) 466 462 DO jl = 1, jpl 467 463 DO jk = 1, nlay_s 468 e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * a rea(:,:) * a_i(:,:,jl) * ht_s(:,:,jl) / ( unit_fac * REAL( nlay_s ))464 e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * a_i(:,:,jl) * ht_s(:,:,jl) / REAL( nlay_s ) 469 465 END DO 470 466 END DO … … 489 485 CALL prt_ctl_info(' - Cell values : ') 490 486 CALL prt_ctl_info(' ~~~~~~~~~~~~~ ') 491 CALL prt_ctl(tab2d_1= area, clinfo1=' lim_thd : cell area :')487 CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_thd : cell area :') 492 488 CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_thd : at_i :') 493 489 CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_thd : vt_i :') … … 520 516 CALL wrk_dealloc( jpi, jpj, zqsr, zqns ) 521 517 518 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 522 519 !------------------------------------------------------------------------------| 523 520 ! 1) Transport of ice between thickness categories. | 524 521 !------------------------------------------------------------------------------| 525 ! Given thermodynamic growth rates, transport ice between 526 ! thickness categories. 522 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 523 524 ! Given thermodynamic growth rates, transport ice between thickness categories. 527 525 IF( jpl > 1 ) CALL lim_itd_th_rem( 1, jpl, kt ) 528 526 ! … … 530 528 CALL lim_var_agg(1) 531 529 530 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 531 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 532 532 !------------------------------------------------------------------------------| 533 533 ! 3) Add frazil ice growing in leads. … … 536 536 CALL lim_var_glo2eqv ! only for info 537 537 538 ! conservation test 539 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 540 538 541 IF(ln_ctl) THEN ! Control print 539 542 CALL prt_ctl_info(' ') 540 543 CALL prt_ctl_info(' - Cell values : ') 541 544 CALL prt_ctl_info(' ~~~~~~~~~~~~~ ') 542 CALL prt_ctl(tab2d_1= area, clinfo1=' lim_itd_th : cell area :')545 CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_itd_th : cell area :') 543 546 CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_th : at_i :') 544 547 CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_th : vt_i :') … … 568 571 ENDIF 569 572 ! 570 ! conservation test571 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)572 573 573 IF( nn_timing == 1 ) CALL timing_stop('limthd') 574 574 … … 614 614 !! ( dA = A/2h dh ) 615 615 !!----------------------------------------------------------------------- 616 INTEGER, INTENT(in) :: kideb, kiut ! bounds for the spatial loop 617 INTEGER :: ji ! dummy loop indices 618 REAL(wp) :: zhi_bef ! ice thickness before thermo 619 REAL(wp) :: zdh_mel ! net melting 616 INTEGER, INTENT(in) :: kideb, kiut ! bounds for the spatial loop 617 INTEGER :: ji ! dummy loop indices 618 REAL(wp) :: zhi_bef ! ice thickness before thermo 619 REAL(wp) :: zdh_mel, zda_mel ! net melting 620 REAL(wp) :: zv ! ice volume 620 621 621 622 DO ji = kideb, kiut 622 zhi_bef = ht_i_1d(ji) - ( dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) ) 623 zdh_mel = MIN( 0._wp, dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) ) 624 IF( zdh_mel < 0._wp ) & 625 & a_i_1d(ji) = MAX( 0._wp, a_i_1d(ji) + a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi10 ) ) ) 623 zdh_mel = MIN( 0._wp, dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) ) 624 IF( zdh_mel < 0._wp ) THEN 625 zv = a_i_1d(ji) * ht_i_1d(ji) 626 ! lateral melting = concentration change 627 zhi_bef = ht_i_1d(ji) - zdh_mel 628 zda_mel = a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi10 ) ) 629 a_i_1d(ji) = MAX( 0._wp, a_i_1d(ji) + zda_mel ) 630 ! adjust thickness 631 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - a_i_1d(ji) + epsi20 ) ) 632 ht_i_1d(ji) = rswitch * zv / MAX( a_i_1d(ji), epsi20 ) 633 ! adjust thickness 634 at_i_1d(ji) = a_i_1d(ji) 635 END IF 626 636 END DO 627 at_i_1d(:) = a_i_1d(:)628 637 629 638 END SUBROUTINE lim_thd_lam … … 665 674 IF ( ( jpl > 1 ) .AND. ( nn_monocat == 1 ) ) THEN 666 675 nn_monocat = 0 667 WRITE(numout, *) 'nn_monocat must be 0 in multi-category case '676 IF(lwp) WRITE(numout, *) ' nn_monocat must be 0 in multi-category case ' 668 677 ENDIF 669 678
Note: See TracChangeset
for help on using the changeset viewer.
