Changeset 2528 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
- Timestamp:
- 2010-12-27T18:33:53+01:00 (13 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90 (modified) (33 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r2477 r2528 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 ice 17 !! lim_thd : thermodynamic of sea ice 18 !! lim_thd_init : initialisation of sea-ice thermodynamic 17 19 !!---------------------------------------------------------------------- 18 20 USE phycst ! physical constants 19 21 USE dom_oce ! ocean space and time domain variables 20 USE ice ! LIM sea-ice variables21 USE par_ice 22 USE ice ! LIM: sea-ice variables 23 USE par_ice ! LIM: sea-ice parameters 22 24 USE sbc_oce ! Surface boundary condition: ocean fields 23 25 USE sbc_ice ! Surface boundary condition: ice fields 24 26 USE thd_ice ! LIM thermodynamic sea-ice variables 25 27 USE dom_ice ! LIM sea-ice domain 26 USE domvvl 27 USE limthd_dif 28 USE limthd_dh 29 USE limthd_sal 30 USE limthd_ent 31 USE limtab 32 USE limvar 28 USE domvvl ! domain: variable volume level 29 USE limthd_dif ! LIM: thermodynamics, vertical diffusion 30 USE limthd_dh ! LIM: thermodynamics, ice and snow thickness variation 31 USE limthd_sal ! LIM: thermodynamics, ice salinity 32 USE limthd_ent ! LIM: thermodynamics, ice enthalpy redistribution 33 USE limtab ! LIM: 1D <==> 2D transformation 34 USE limvar ! LIM: sea-ice variables 35 USE lbclnk ! lateral boundary condition - MPP links 36 USE lib_mpp ! MPP library 33 37 USE in_out_manager ! I/O manager 34 38 USE prtctl ! Print control 35 USE lbclnk36 USE lib_mpp37 39 38 40 IMPLICIT NONE 39 41 PRIVATE 40 42 41 PUBLIC lim_thd ! called by lim_step 42 43 REAL(wp) :: epsi20 = 1e-20 ! constant values 44 REAL(wp) :: epsi16 = 1e-16 ! 45 REAL(wp) :: epsi06 = 1e-06 ! 46 REAL(wp) :: epsi04 = 1e-04 ! 47 REAL(wp) :: zzero = 0.e0 ! 48 REAL(wp) :: zone = 1.e0 ! 43 PUBLIC lim_thd ! called by limstp module 44 PUBLIC lim_thd_init ! called by iceini module 45 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 ! 49 52 50 53 !! * Substitutions … … 52 55 # include "vectopt_loop_substitute.h90" 53 56 !!---------------------------------------------------------------------- 54 !! NEMO/LIM 3.2 , UCL-LOCEAN-IPSL (2009)57 !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010) 55 58 !! $Id$ 56 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)59 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 57 60 !!---------------------------------------------------------------------- 58 59 61 CONTAINS 60 62 … … 84 86 REAL(wp) :: zfric_umax = 2e-02 ! upper bound for the friction velocity 85 87 REAL(wp) :: zinda, zindb, zthsnice, zfric_u ! temporary scalar 86 REAL(wp) :: zfntlat, zpareff 88 REAL(wp) :: zfntlat, zpareff ! - - 87 89 REAL(wp) :: zeps, zareamin, zcoef 88 90 REAL(wp), DIMENSION(jpi,jpj) :: zqlbsbq ! link with lead energy budget qldif 89 91 !!------------------------------------------------------------------- 90 91 92 92 93 !------------------------------------------------------------------------------! … … 99 100 !-------------------- 100 101 ! Change the units of heat content; from global units to J.m3 101 102 102 DO jl = 1, jpl 103 103 DO jk = 1, nlay_i … … 166 166 pfrld(ji,jj) = 1.0 - at_i(ji,jj) 167 167 zinda = 1.0 - MAX( zzero , SIGN( zone , - ( 1.0 - pfrld(ji,jj) ) ) ) 168 168 ! 169 169 ! ! solar irradiance transmission at the mixed layer bottom and used in the lead heat budget 170 170 ! ! practically no "direct lateral ablation" … … 181 181 qdtcn(ji,jj) = zindb * fdtcn(ji,jj) * (1.0 - at_i(ji,jj)) * rdt_ice 182 182 ! 183 184 ! still need to be updated : fdtcn !!!! 185 ! !-- Lead heat budget (part 1, next one is in limthd_dh 186 ! !-- qldif -- (or qldif_1d in 1d routines) 187 qldif(ji,jj) = tms(ji,jj) * rdt_ice * ( & 188 & pfrld(ji,jj) * ( qsr(ji,jj) & ! solar heat 189 & + qns(ji,jj) & ! non solar heat 190 & + fdtcn(ji,jj) & ! turbulent ice-ocean heat 191 & + fsbbq(ji,jj) * ( 1.0 - zindb ) ) & ! residual heat from previous step 183 ! !-- Lead heat budget, qldif (part 1, next one is in limthd_dh) 184 ! ! caution: exponent betas used as more snow can fallinto leads 185 qldif(ji,jj) = tms(ji,jj) * rdt_ice * ( & 186 & pfrld(ji,jj) * ( qsr(ji,jj) & ! solar heat 187 & + qns(ji,jj) & ! non solar heat 188 & + fdtcn(ji,jj) & ! turbulent ice-ocean heat 189 & + fsbbq(ji,jj) * ( 1.0 - zindb ) ) & ! residual heat from previous step 192 190 & - pfrld(ji,jj)**betas * sprecip(ji,jj) * lfus ) ! latent heat of sprecip melting 193 191 ! 194 192 ! Positive heat budget is used for bottom ablation 195 193 zfntlat = 1.0 - MAX( zzero , SIGN( zone , - qldif(ji,jj) ) ) … … 198 196 != 0 if ice and positive heat budget and 1 if one of those two is false 199 197 zqlbsbq(ji,jj) = qldif(ji,jj) * ( 1.0 - zpareff ) / MAX( at_i(ji,jj) * rdt_ice , epsi16 ) 200 198 ! 201 199 ! Heat budget of the lead, energy transferred from ice to ocean 202 200 qldif (ji,jj) = zpareff * qldif(ji,jj) 203 201 qdtcn (ji,jj) = zpareff * qdtcn(ji,jj) 204 202 ! 205 203 ! Energy needed to bring ocean surface layer until its freezing (qcmif, limflx) 206 204 qcmif (ji,jj) = rau0 * rcp * fse3t(ji,jj,1) * ( t_bo(ji,jj) - (sst_m(ji,jj) + rt0) ) * ( 1. - zinda ) 207 205 ! 208 206 ! oceanic heat flux (limthd_dh) 209 207 fbif (ji,jj) = zindb * ( fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi20 ) + fdtcn(ji,jj) ) … … 223 221 ENDIF 224 222 225 zareamin = 1. 0e-10223 zareamin = 1.e-10 226 224 nbpb = 0 227 225 DO jj = 1, jpj … … 359 357 CALL tab_1d_2d( nbpb, rdvonif, npb, dvnbq_1d (1:nbpb), jpi, jpj ) 360 358 CALL tab_1d_2d( nbpb, fseqv , npb, fseqv_1d (1:nbpb), jpi, jpj ) 361 359 ! 362 360 IF( num_sal == 2 ) THEN 363 361 CALL tab_1d_2d( nbpb, fsbri, npb, fsbri_1d(1:nbpb), jpi, jpj ) 364 362 CALL tab_1d_2d( nbpb, fhbri, npb, fhbri_1d(1:nbpb), jpi, jpj ) 365 363 ENDIF 366 364 ! 367 365 !+++++ 368 366 !temporary stuff for a dummy version … … 375 373 CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qnsr_ice_1d(1:nbpb), jpi, jpj) 376 374 !+++++ 377 375 ! 378 376 IF( lk_mpp ) CALL mpp_comm_free( ncomm_ice ) !RB necessary ?? 379 377 ENDIF … … 388 386 ! 5.1) Ice heat content 389 387 !------------------------ 390 391 388 ! Enthalpies are global variables we have to readjust the units 392 389 zcoef = 1.e0 / ( unit_fac * REAL(nlay_i) ) … … 401 398 ! 5.2) Snow heat content 402 399 !------------------------ 403 404 400 ! Enthalpies are global variables we have to readjust the units 405 401 zcoef = 1.e0 / ( unit_fac * REAL(nlay_s) ) … … 424 420 IF( con_i ) fbif(:,:) = fbif(:,:) + zqlbsbq(:,:) 425 421 426 IF(ln_ctl) THEN ! Control print422 IF(ln_ctl) THEN ! Control print 427 423 CALL prt_ctl_info(' ') 428 424 CALL prt_ctl_info(' - Cell values : ') … … 454 450 END DO 455 451 END DO 456 457 452 ENDIF 458 453 ! 459 454 END SUBROUTINE lim_thd 460 455 … … 524 519 END DO 525 520 ! layer by layer 526 dq_i_layer(:,:) 521 dq_i_layer(:,:) = q_i_layer_fin(:,:) - q_i_layer_in(:,:) 527 522 528 523 !---------------------------------------- 529 524 ! Atmospheric heat flux, ice heat budget 530 525 !---------------------------------------- 531 532 DO ji = kideb, kiut 533 zji = MOD( npb(ji) - 1, jpi ) + 1 534 zjj = ( npb(ji) - 1 ) / jpi + 1 535 536 fatm(ji,jl) = qnsr_ice_1d(ji) + (1.0-i0(ji))*qsr_ice_1d(ji) 537 538 sum_fluxq(ji,jl) = fc_su(ji) - fc_bo_i(ji) + qsr_ice_1d(ji)*i0(ji) - fstroc(zji,zjj,jl) 526 DO ji = kideb, kiut 527 zji = MOD( npb(ji) - 1 , jpi ) + 1 528 zjj = ( npb(ji) - 1 ) / jpi + 1 529 fatm (ji,jl) = qnsr_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) 530 sum_fluxq(ji,jl) = fc_su(ji) - fc_bo_i(ji) + qsr_ice_1d(ji) * i0(ji) - fstroc(zji,zjj,jl) 539 531 END DO 540 532 … … 542 534 ! Conservation error 543 535 !-------------------- 544 545 536 DO ji = kideb, kiut 546 537 cons_error(ji,jl) = ABS( dq_i(ji,jl) / rdt_ice + sum_fluxq(ji,jl) ) 547 538 END DO 548 539 549 numce = 0540 numce = 0 550 541 meance = 0.0 551 542 DO ji = kideb, kiut … … 554 545 meance = meance + cons_error(ji,jl) 555 546 ENDIF 556 END DO557 IF (numce .GT. 0 ) meance = meance / numce547 END DO 548 IF( numce .GT. 0 ) meance = meance / numce 558 549 559 550 WRITE(numout,*) ' Maximum tolerated conservation error : ', max_cons_err … … 565 556 ! Surface error due to imbalance between Fatm and Fcsu 566 557 !------------------------------------------------------- 567 numce = 0 .0558 numce = 0 568 559 meance = 0.0 569 560 570 561 DO ji = kideb, kiut 571 562 surf_error(ji,jl) = ABS ( fatm(ji,jl) - fc_su(ji) ) 572 IF ( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. & 573 max_surf_err ) ) THEN 563 IF( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. max_surf_err ) ) THEN 574 564 numce = numce + 1 575 565 meance = meance + surf_error(ji,jl) 576 566 ENDIF 577 567 ENDDO 578 IF (numce .GT. 0 ) meance = meance / numce568 IF( numce .GT. 0 ) meance = meance / numce 579 569 580 570 WRITE(numout,*) ' Maximum tolerated surface error : ', max_surf_err … … 590 580 ! Write ice state in case of big errors 591 581 !--------------------------------------- 592 593 582 DO ji = kideb, kiut 594 583 IF ( ( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. max_surf_err ) ) .OR. & … … 596 585 zji = MOD( npb(ji) - 1, jpi ) + 1 597 586 zjj = ( npb(ji) - 1 ) / jpi + 1 598 587 ! 599 588 WRITE(numout,*) ' alerte 1 ' 600 589 WRITE(numout,*) ' Untolerated conservation / surface error after ' … … 612 601 ! WRITE(numout,*) ' qt_i_fin : ', qt_i_fin(ji,jl) 613 602 ! WRITE(numout,*) ' qt_s_fin : ', qt_s_fin(ji,jl) 614 ! WRITE(numout,*) ' qt : ', qt_i_fin(ji,jl) + & 615 ! qt_s_fin(ji,jl) 603 ! WRITE(numout,*) ' qt : ', qt_i_fin(ji,jl) + qt_s_fin(ji,jl) 616 604 WRITE(numout,*) ' ht_i : ', ht_i_b(ji) 617 605 WRITE(numout,*) ' ht_s : ', ht_s_b(ji) … … 640 628 WRITE(numout,*) 641 629 WRITE(numout,*) ' Layer by layer ... ' 642 WRITE(numout,*) ' dq_snow : ', ( qt_s_fin(ji,jl) - & 643 qt_s_in(ji,jl) ) & 644 / rdt_ice 645 WRITE(numout,*) ' dfc_snow : ', fc_s(ji,1) - & 646 fc_s(ji,0) 630 WRITE(numout,*) ' dq_snow : ', ( qt_s_fin(ji,jl) - qt_s_in(ji,jl) ) / rdt_ice 631 WRITE(numout,*) ' dfc_snow : ', fc_s(ji,1) - fc_s(ji,0) 647 632 DO jk = 1, nlay_i 648 633 WRITE(numout,*) ' layer : ', jk 649 634 WRITE(numout,*) ' dq_ice : ', dq_i_layer(ji,jk) / rdt_ice 650 635 WRITE(numout,*) ' radab : ', radab(ji,jk) 651 WRITE(numout,*) ' dfc_i : ', fc_i(ji,jk) - & 652 fc_i(ji,jk-1) 653 WRITE(numout,*) ' tot f : ', fc_i(ji,jk) - & 654 fc_i(ji,jk-1) - radab(ji,jk) 636 WRITE(numout,*) ' dfc_i : ', fc_i(ji,jk) - fc_i(ji,jk-1) 637 WRITE(numout,*) ' tot f : ', fc_i(ji,jk) - fc_i(ji,jk-1) - radab(ji,jk) 655 638 END DO 656 639 … … 662 645 663 646 664 SUBROUTINE lim_thd_con_dh( kideb,kiut,jl)647 SUBROUTINE lim_thd_con_dh( kideb, kiut, jl ) 665 648 !!----------------------------------------------------------------------- 666 649 !! *** ROUTINE lim_thd_con_dh *** 667 650 !! 668 651 !! ** Purpose : Test energy conservation after enthalpy redistr. 669 !!670 !! history :671 !! 9.9 ! 07-04 (M.Vancoppenolle) original code672 652 !!----------------------------------------------------------------------- 673 653 INTEGER, INTENT(in) :: & … … 745 725 ! Write ice state in case of big errors 746 726 !--------------------------------------- 747 748 727 DO ji = kideb, kiut 749 728 IF ( cons_error(ji,jl) .GT. max_cons_err ) THEN 750 729 zji = MOD( npb(ji) - 1, jpi ) + 1 751 730 zjj = ( npb(ji) - 1 ) / jpi + 1 752 731 ! 753 732 WRITE(numout,*) ' alerte 1 - category : ', jl 754 733 WRITE(numout,*) ' Untolerated conservation error after limthd_ent ' … … 771 750 WRITE(numout,*) ' qt_s_in : ', qt_s_in(ji,jl) / rdt_ice 772 751 WRITE(numout,*) ' qt_i_in : ', qt_i_in(ji,jl) / rdt_ice 773 WRITE(numout,*) ' qt_in : ', ( qt_i_in(ji,jl) + & 774 qt_s_in(ji,jl) ) / rdt_ice 752 WRITE(numout,*) ' qt_in : ', ( qt_i_in(ji,jl) + qt_s_in(ji,jl) ) / rdt_ice 775 753 WRITE(numout,*) ' qt_s_fin : ', qt_s_fin(ji,jl) / rdt_ice 776 754 WRITE(numout,*) ' qt_i_fin : ', qt_i_fin(ji,jl) / rdt_ice 777 WRITE(numout,*) ' qt_fin : ', ( qt_i_fin(ji,jl) + & 778 qt_s_fin(ji,jl) ) / rdt_ice 755 WRITE(numout,*) ' qt_fin : ', ( qt_i_fin(ji,jl) + qt_s_fin(ji,jl) ) / rdt_ice 779 756 WRITE(numout,*) ' * ' 780 757 WRITE(numout,*) ' Ice variables --- : ' … … 785 762 WRITE(numout,*) ' dh_i_surf : ', dh_i_surf(ji) 786 763 WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji) 787 788 764 ENDIF 789 765 ! … … 800 776 !! 801 777 !! ** Method : Formula (Bitz and Lipscomb, 1999) 802 !!803 778 !!------------------------------------------------------------------- 804 779 INTEGER, INTENT(in) :: kideb, kiut ! bounds for the spatial loop … … 827 802 828 803 SUBROUTINE lim_thd_init 829 830 804 !!----------------------------------------------------------------------- 831 805 !! *** ROUTINE lim_thd_init *** … … 838 812 !! 839 813 !! ** input : Namelist namicether 840 814 !!------------------------------------------------------------------- 841 815 NAMELIST/namicethd/ hmelt , hiccrit, fraz_swi, maxfrazb, vfrazb, Cfrazb, & 842 816 & hicmin, hiclim, amax , & … … 845 819 & kappa_i, nconv_i_thd, maxer_i_thd, thcon_i_swi 846 820 !!------------------------------------------------------------------- 847 821 ! 848 822 IF(lwp) THEN 849 823 WRITE(numout,*) … … 851 825 WRITE(numout,*) '~~~~~~~' 852 826 ENDIF 853 827 ! 854 828 REWIND( numnam_ice ) ! read Namelist numnam_ice 855 829 READ ( numnam_ice , namicethd ) 856 830 ! 857 831 IF(lwp) THEN ! control print 858 832 WRITE(numout,*) … … 892 866 #else 893 867 !!---------------------------------------------------------------------- 894 !! Default option 868 !! Default option Dummy module NO LIM3 sea-ice model 895 869 !!---------------------------------------------------------------------- 896 CONTAINS897 SUBROUTINE lim_thd ! Empty routine898 END SUBROUTINE lim_thd899 SUBROUTINE lim_thd_con_dif900 END SUBROUTINE lim_thd_con_dif901 870 #endif 902 871
Note: See TracChangeset
for help on using the changeset viewer.
