- Timestamp:
- 2010-11-13T16:01:26+01:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r2287 r2383 7 7 !! 2.0 ! 2002-07 (C. Ethe, G. Madec) LIM-2 (F90 rewriting) 8 8 !! 3.0 ! 2005-11 (M. Vancoppenolle) LIM-3 : Multi-layer thermodynamics + salinity variations 9 !! - ! 2007-04 (M. Vancoppenolle) add lim_thd_glohec and lim_thd_con_dif9 !! - ! 2007-04 (M. Vancoppenolle) add lim_thd_glohec, lim_thd_con_dh and lim_thd_con_dif 10 10 !! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in rdmsnif 11 !! 3.3 ! 2010-11 (G. Madec) corrected snow melting heat (due to factor betas) 11 12 !!---------------------------------------------------------------------- 12 13 #if defined key_lim3 … … 14 15 !! 'key_lim3' LIM3 sea-ice model 15 16 !!---------------------------------------------------------------------- 16 !! lim_thd : thermodynamic of sea ice17 !! lim_thd_init : initialisation of sea-ice thermodynamic17 !! lim_thd : thermodynamic of sea ice 18 !! lim_thd_init : initialisation of sea-ice thermodynamic 18 19 !!---------------------------------------------------------------------- 19 20 USE phycst ! physical constants 20 21 USE dom_oce ! ocean space and time domain variables 21 USE ice ! LIM sea-ice variables22 USE par_ice 22 USE ice ! LIM: sea-ice variables 23 USE par_ice ! LIM: sea-ice parameters 23 24 USE sbc_oce ! Surface boundary condition: ocean fields 24 25 USE sbc_ice ! Surface boundary condition: ice fields 25 26 USE thd_ice ! LIM thermodynamic sea-ice variables 26 27 USE dom_ice ! LIM sea-ice domain 27 USE domvvl 28 USE iceini 29 USE limthd_dif 30 USE limthd_dh 31 USE limthd_sal 32 USE limthd_ent 33 USE limtab 34 USE limvar 28 USE domvvl ! domain: variable volume level 29 USE iceini ! LIM: sea-ice initialization 30 USE limthd_dif ! LIM: thermodynamics, vertical diffusion 31 USE limthd_dh ! LIM: thermodynamics, ice and snow thickness variation 32 USE limthd_sal ! LIM: thermodynamics, ice salinity 33 USE limthd_ent ! LIM: thermodynamics, ice enthalpy redistribution 34 USE limtab ! LIM: 1D <==> 2D transformation 35 USE limvar ! LIM: ??? 36 USE lbclnk ! lateral boundary condition - MPP links 37 USE lib_mpp ! MPP library 35 38 USE in_out_manager ! I/O manager 36 39 USE prtctl ! Print control 37 USE lbclnk38 USE lib_mpp39 40 40 41 IMPLICIT NONE … … 43 44 PUBLIC lim_thd ! called by lim_step 44 45 45 REAL(wp) :: epsi20 = 1e-20 ! constant values46 REAL(wp) :: epsi16 = 1e-16 !47 REAL(wp) :: epsi06 = 1e-06 !48 REAL(wp) :: epsi04 = 1e-04 !49 REAL(wp) :: zzero = 0. e0!50 REAL(wp) :: zone = 1. e0!46 REAL(wp) :: epsi20 = 1e-20_wp ! constant values 47 REAL(wp) :: epsi16 = 1e-16_wp ! 48 REAL(wp) :: epsi06 = 1e-06_wp ! 49 REAL(wp) :: epsi04 = 1e-04_wp ! 50 REAL(wp) :: zzero = 0._wp ! 51 REAL(wp) :: zone = 1._wp ! 51 52 52 53 !! * Substitutions … … 56 57 !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010) 57 58 !! $Id$ 58 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)59 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 59 60 !!---------------------------------------------------------------------- 60 61 61 CONTAINS 62 62 … … 86 86 REAL(wp) :: zfric_umax = 2e-02 ! upper bound for the friction velocity 87 87 REAL(wp) :: zinda, zindb, zthsnice, zfric_u ! temporary scalar 88 REAL(wp) :: zfn sol, zfontn, zfntlat, zpareff ! - -88 REAL(wp) :: zfntlat, zpareff ! - - 89 89 REAL(wp) :: zeps, zareamin, zcoef 90 90 REAL(wp), DIMENSION(jpi,jpj) :: zqlbsbq ! link with lead energy budget qldif … … 104 104 !-------------------- 105 105 ! Change the units of heat content; from global units to J.m3 106 107 106 DO jl = 1, jpl 108 107 DO jk = 1, nlay_i … … 171 170 pfrld(ji,jj) = 1.0 - at_i(ji,jj) 172 171 zinda = 1.0 - MAX( zzero , SIGN( zone , - ( 1.0 - pfrld(ji,jj) ) ) ) 173 172 ! 174 173 ! ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget 175 174 ! ! practically no "direct lateral ablation" … … 186 185 qdtcn(ji,jj) = zindb * fdtcn(ji,jj) * (1.0 - at_i(ji,jj)) * rdt_ice 187 186 ! 188 189 ! still need to be updated : fdtcn !!!! 190 ! !-- Lead heat budget (part 1, next one is in limthd_dh 191 ! !-- qldif -- (or qldif_1d in 1d routines) 192 zfontn = sprecip(ji,jj) * lfus ! energy of melting 193 zfnsol = qns(ji,jj) ! total non solar flux 194 qldif(ji,jj) = tms(ji,jj) * ( qsr(ji,jj) & 195 & + zfnsol + fdtcn(ji,jj) - zfontn & 196 & + ( 1.0 - zindb ) * fsbbq(ji,jj) ) & 197 & * ( 1.0 - at_i(ji,jj) ) * rdt_ice 198 187 ! !-- Lead heat budget, qldif (part 1, next one is in limthd_dh) 188 ! ! caution: exponent betas used as more snow can fallinto leads 189 qldif(ji,jj) = tms(ji,jj) * rdt_ice * ( & 190 & pfrld(ji,jj) * ( qsr(ji,jj) & ! solar heat 191 & + qns(ji,jj) & ! non solar heat 192 & + fdtcn(ji,jj) & ! turbulent ice-ocean heat 193 & + fsbbq(ji,jj) * ( 1.0 - zindb ) ) & ! residual heat from previous step 194 & - pfrld(ji,jj)**betas * sprecip(ji,jj) * lfus ) ! latent heat of sprecip melting 195 ! 199 196 ! Positive heat budget is used for bottom ablation 200 197 zfntlat = 1.0 - MAX( zzero , SIGN( zone , - qldif(ji,jj) ) ) … … 203 200 != 0 if ice and positive heat budget and 1 if one of those two is false 204 201 zqlbsbq(ji,jj) = qldif(ji,jj) * ( 1.0 - zpareff ) / MAX( at_i(ji,jj) * rdt_ice , epsi16 ) 205 202 ! 206 203 ! Heat budget of the lead, energy transferred from ice to ocean 207 204 qldif (ji,jj) = zpareff * qldif(ji,jj) 208 205 qdtcn (ji,jj) = zpareff * qdtcn(ji,jj) 209 206 ! 210 207 ! Energy needed to bring ocean surface layer until its freezing (qcmif, limflx) 211 208 qcmif (ji,jj) = rau0 * rcp * fse3t(ji,jj,1) * ( t_bo(ji,jj) - (sst_m(ji,jj) + rt0) ) * ( 1. - zinda ) 212 209 ! 213 210 ! oceanic heat flux (limthd_dh) 214 211 fbif (ji,jj) = zindb * ( fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi20 ) + fdtcn(ji,jj) ) … … 228 225 ENDIF 229 226 230 zareamin = 1. 0e-10227 zareamin = 1.e-10 231 228 nbpb = 0 232 229 DO jj = 1, jpj … … 364 361 CALL tab_1d_2d( nbpb, rdvonif, npb, dvnbq_1d (1:nbpb), jpi, jpj ) 365 362 CALL tab_1d_2d( nbpb, fseqv , npb, fseqv_1d (1:nbpb), jpi, jpj ) 366 363 ! 367 364 IF( num_sal == 2 ) THEN 368 365 CALL tab_1d_2d( nbpb, fsbri, npb, fsbri_1d(1:nbpb), jpi, jpj ) 369 366 CALL tab_1d_2d( nbpb, fhbri, npb, fhbri_1d(1:nbpb), jpi, jpj ) 370 367 ENDIF 371 368 ! 372 369 !+++++ 373 370 !temporary stuff for a dummy version … … 380 377 CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qnsr_ice_1d(1:nbpb), jpi, jpj) 381 378 !+++++ 382 379 ! 383 380 IF( lk_mpp ) CALL mpp_comm_free( ncomm_ice ) !RB necessary ?? 384 381 ENDIF … … 393 390 ! 5.1) Ice heat content 394 391 !------------------------ 395 396 392 ! Enthalpies are global variables we have to readjust the units 397 393 zcoef = 1.e0 / ( unit_fac * REAL(nlay_i) ) … … 406 402 ! 5.2) Snow heat content 407 403 !------------------------ 408 409 404 ! Enthalpies are global variables we have to readjust the units 410 405 zcoef = 1.e0 / ( unit_fac * REAL(nlay_s) ) … … 429 424 IF( con_i ) fbif(:,:) = fbif(:,:) + zqlbsbq(:,:) 430 425 431 IF(ln_ctl) THEN ! Control print426 IF(ln_ctl) THEN ! Control print 432 427 CALL prt_ctl_info(' ') 433 428 CALL prt_ctl_info(' - Cell values : ') … … 459 454 END DO 460 455 END DO 461 462 456 ENDIF 463 457 ! 464 458 END SUBROUTINE lim_thd 465 459 … … 529 523 END DO 530 524 ! layer by layer 531 dq_i_layer(:,:) 525 dq_i_layer(:,:) = q_i_layer_fin(:,:) - q_i_layer_in(:,:) 532 526 533 527 !---------------------------------------- 534 528 ! Atmospheric heat flux, ice heat budget 535 529 !---------------------------------------- 536 537 DO ji = kideb, kiut 538 zji = MOD( npb(ji) - 1, jpi ) + 1 539 zjj = ( npb(ji) - 1 ) / jpi + 1 540 541 fatm(ji,jl) = qnsr_ice_1d(ji) + (1.0-i0(ji))*qsr_ice_1d(ji) 542 543 sum_fluxq(ji,jl) = fc_su(ji) - fc_bo_i(ji) + qsr_ice_1d(ji)*i0(ji) - fstroc(zji,zjj,jl) 530 DO ji = kideb, kiut 531 zji = MOD( npb(ji) - 1 , jpi ) + 1 532 zjj = ( npb(ji) - 1 ) / jpi + 1 533 fatm (ji,jl) = qnsr_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) 534 sum_fluxq(ji,jl) = fc_su(ji) - fc_bo_i(ji) + qsr_ice_1d(ji) * i0(ji) - fstroc(zji,zjj,jl) 544 535 END DO 545 536 … … 547 538 ! Conservation error 548 539 !-------------------- 549 550 540 DO ji = kideb, kiut 551 541 cons_error(ji,jl) = ABS( dq_i(ji,jl) / rdt_ice + sum_fluxq(ji,jl) ) 552 542 END DO 553 543 554 numce = 0544 numce = 0 555 545 meance = 0.0 556 546 DO ji = kideb, kiut … … 559 549 meance = meance + cons_error(ji,jl) 560 550 ENDIF 561 END DO562 IF (numce .GT. 0 ) meance = meance / numce551 END DO 552 IF( numce .GT. 0 ) meance = meance / numce 563 553 564 554 WRITE(numout,*) ' Maximum tolerated conservation error : ', max_cons_err … … 570 560 ! Surface error due to imbalance between Fatm and Fcsu 571 561 !------------------------------------------------------- 572 numce = 0 .0562 numce = 0 573 563 meance = 0.0 574 564 575 565 DO ji = kideb, kiut 576 566 surf_error(ji,jl) = ABS ( fatm(ji,jl) - fc_su(ji) ) 577 IF ( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. & 578 max_surf_err ) ) THEN 567 IF( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. max_surf_err ) ) THEN 579 568 numce = numce + 1 580 569 meance = meance + surf_error(ji,jl) 581 570 ENDIF 582 571 ENDDO 583 IF (numce .GT. 0 ) meance = meance / numce572 IF( numce .GT. 0 ) meance = meance / numce 584 573 585 574 WRITE(numout,*) ' Maximum tolerated surface error : ', max_surf_err … … 595 584 ! Write ice state in case of big errors 596 585 !--------------------------------------- 597 598 586 DO ji = kideb, kiut 599 587 IF ( ( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. max_surf_err ) ) .OR. & … … 601 589 zji = MOD( npb(ji) - 1, jpi ) + 1 602 590 zjj = ( npb(ji) - 1 ) / jpi + 1 603 591 ! 604 592 WRITE(numout,*) ' alerte 1 ' 605 593 WRITE(numout,*) ' Untolerated conservation / surface error after ' … … 617 605 ! WRITE(numout,*) ' qt_i_fin : ', qt_i_fin(ji,jl) 618 606 ! WRITE(numout,*) ' qt_s_fin : ', qt_s_fin(ji,jl) 619 ! WRITE(numout,*) ' qt : ', qt_i_fin(ji,jl) + & 620 ! qt_s_fin(ji,jl) 607 ! WRITE(numout,*) ' qt : ', qt_i_fin(ji,jl) + qt_s_fin(ji,jl) 621 608 WRITE(numout,*) ' ht_i : ', ht_i_b(ji) 622 609 WRITE(numout,*) ' ht_s : ', ht_s_b(ji) … … 645 632 WRITE(numout,*) 646 633 WRITE(numout,*) ' Layer by layer ... ' 647 WRITE(numout,*) ' dq_snow : ', ( qt_s_fin(ji,jl) - & 648 qt_s_in(ji,jl) ) & 649 / rdt_ice 650 WRITE(numout,*) ' dfc_snow : ', fc_s(ji,1) - & 651 fc_s(ji,0) 634 WRITE(numout,*) ' dq_snow : ', ( qt_s_fin(ji,jl) - qt_s_in(ji,jl) ) / rdt_ice 635 WRITE(numout,*) ' dfc_snow : ', fc_s(ji,1) - fc_s(ji,0) 652 636 DO jk = 1, nlay_i 653 637 WRITE(numout,*) ' layer : ', jk 654 638 WRITE(numout,*) ' dq_ice : ', dq_i_layer(ji,jk) / rdt_ice 655 639 WRITE(numout,*) ' radab : ', radab(ji,jk) 656 WRITE(numout,*) ' dfc_i : ', fc_i(ji,jk) - & 657 fc_i(ji,jk-1) 658 WRITE(numout,*) ' tot f : ', fc_i(ji,jk) - & 659 fc_i(ji,jk-1) - radab(ji,jk) 640 WRITE(numout,*) ' dfc_i : ', fc_i(ji,jk) - fc_i(ji,jk-1) 641 WRITE(numout,*) ' tot f : ', fc_i(ji,jk) - fc_i(ji,jk-1) - radab(ji,jk) 660 642 END DO 661 643 … … 667 649 668 650 669 SUBROUTINE lim_thd_con_dh( kideb,kiut,jl)651 SUBROUTINE lim_thd_con_dh( kideb, kiut, jl ) 670 652 !!----------------------------------------------------------------------- 671 653 !! *** ROUTINE lim_thd_con_dh *** 672 654 !! 673 655 !! ** Purpose : Test energy conservation after enthalpy redistr. 674 !!675 !! history :676 !! 9.9 ! 07-04 (M.Vancoppenolle) original code677 656 !!----------------------------------------------------------------------- 678 657 INTEGER, INTENT(in) :: & … … 750 729 ! Write ice state in case of big errors 751 730 !--------------------------------------- 752 753 731 DO ji = kideb, kiut 754 732 IF ( cons_error(ji,jl) .GT. max_cons_err ) THEN 755 733 zji = MOD( npb(ji) - 1, jpi ) + 1 756 734 zjj = ( npb(ji) - 1 ) / jpi + 1 757 735 ! 758 736 WRITE(numout,*) ' alerte 1 - category : ', jl 759 737 WRITE(numout,*) ' Untolerated conservation error after limthd_ent ' … … 776 754 WRITE(numout,*) ' qt_s_in : ', qt_s_in(ji,jl) / rdt_ice 777 755 WRITE(numout,*) ' qt_i_in : ', qt_i_in(ji,jl) / rdt_ice 778 WRITE(numout,*) ' qt_in : ', ( qt_i_in(ji,jl) + & 779 qt_s_in(ji,jl) ) / rdt_ice 756 WRITE(numout,*) ' qt_in : ', ( qt_i_in(ji,jl) + qt_s_in(ji,jl) ) / rdt_ice 780 757 WRITE(numout,*) ' qt_s_fin : ', qt_s_fin(ji,jl) / rdt_ice 781 758 WRITE(numout,*) ' qt_i_fin : ', qt_i_fin(ji,jl) / rdt_ice 782 WRITE(numout,*) ' qt_fin : ', ( qt_i_fin(ji,jl) + & 783 qt_s_fin(ji,jl) ) / rdt_ice 759 WRITE(numout,*) ' qt_fin : ', ( qt_i_fin(ji,jl) + qt_s_fin(ji,jl) ) / rdt_ice 784 760 WRITE(numout,*) ' * ' 785 761 WRITE(numout,*) ' Ice variables --- : ' … … 790 766 WRITE(numout,*) ' dh_i_surf : ', dh_i_surf(ji) 791 767 WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji) 792 793 768 ENDIF 794 769 ! … … 805 780 !! 806 781 !! ** Method : Formula (Bitz and Lipscomb, 1999) 807 !!808 782 !!------------------------------------------------------------------- 809 783 INTEGER, INTENT(in) :: kideb, kiut ! bounds for the spatial loop … … 832 806 833 807 SUBROUTINE lim_thd_init 834 835 808 !!----------------------------------------------------------------------- 836 809 !! *** ROUTINE lim_thd_init *** … … 843 816 !! 844 817 !! ** input : Namelist namicether 845 818 !!------------------------------------------------------------------- 846 819 NAMELIST/namicethd/ hmelt , hiccrit, fraz_swi, maxfrazb, vfrazb, Cfrazb, & 847 820 & hicmin, hiclim, amax , &
Note: See TracChangeset
for help on using the changeset viewer.
