Changeset 4688 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
- Timestamp:
- 2014-06-25T01:39:59+02:00 (10 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90 (modified) (13 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r4333 r4688 25 25 USE wrk_nemo ! work arrays 26 26 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 27 USE cpl_oasis3, ONLY : lk_cpl 27 28 28 29 IMPLICIT NONE … … 31 32 PUBLIC lim_thd_dif ! called by lim_thd 32 33 33 REAL(wp) :: epsi10 =1.e-10_wp !34 REAL(wp) :: epsi10 = 1.e-10_wp ! 34 35 !!---------------------------------------------------------------------- 35 36 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) … … 39 40 CONTAINS 40 41 41 SUBROUTINE lim_thd_dif( kideb , kiut , jl)42 SUBROUTINE lim_thd_dif( kideb , kiut ) 42 43 !!------------------------------------------------------------------ 43 44 !! *** ROUTINE lim_thd_dif *** … … 91 92 !! (04-2007) Energy conservation tested by M. Vancoppenolle 92 93 !!------------------------------------------------------------------ 93 INTEGER , INTENT (in) :: kideb ! Start point on which the the computation is applied 94 INTEGER , INTENT (in) :: kiut ! End point on which the the computation is applied 95 INTEGER , INTENT (in) :: jl ! Category number 94 INTEGER , INTENT(in) :: kideb, kiut ! Start/End point on which the the computation is applied 96 95 97 96 !! * Local variables … … 102 101 INTEGER :: nconv ! number of iterations in iterative procedure 103 102 INTEGER :: minnumeqmin, maxnumeqmax 104 INTEGER, DIMENSION(kiut) :: numeqmin ! reference number of top equation105 INTEGER, DIMENSION(kiut) :: numeqmax ! reference number of bottom equation106 INTEGER, DIMENSION(kiut) :: isnow ! switch for presence (1) or absence (0) of snow103 INTEGER, POINTER, DIMENSION(:) :: numeqmin ! reference number of top equation 104 INTEGER, POINTER, DIMENSION(:) :: numeqmax ! reference number of bottom equation 105 INTEGER, POINTER, DIMENSION(:) :: isnow ! switch for presence (1) or absence (0) of snow 107 106 REAL(wp) :: zg1s = 2._wp ! for the tridiagonal system 108 107 REAL(wp) :: zg1 = 2._wp ! … … 111 110 REAL(wp) :: zraext_s = 1.e+8_wp ! extinction coefficient of radiation in the snow 112 111 REAL(wp) :: zkimin = 0.10_wp ! minimum ice thermal conductivity 112 REAL(wp) :: ztsu_err = 1.e-5_wp ! range around which t_su is considered as 0°C 113 113 REAL(wp) :: ztmelt_i ! ice melting temperature 114 114 REAL(wp) :: zerritmax ! current maximal error on temperature 115 REAL(wp), DIMENSION(kiut) :: ztfs ! ice melting point 116 REAL(wp), DIMENSION(kiut) :: ztsuold ! old surface temperature (before the iterative procedure ) 117 REAL(wp), DIMENSION(kiut) :: ztsuoldit ! surface temperature at previous iteration 118 REAL(wp), DIMENSION(kiut) :: zh_i ! ice layer thickness 119 REAL(wp), DIMENSION(kiut) :: zh_s ! snow layer thickness 120 REAL(wp), DIMENSION(kiut) :: zfsw ! solar radiation absorbed at the surface 121 REAL(wp), DIMENSION(kiut) :: zf ! surface flux function 122 REAL(wp), DIMENSION(kiut) :: dzf ! derivative of the surface flux function 123 REAL(wp), DIMENSION(kiut) :: zerrit ! current error on temperature 124 REAL(wp), DIMENSION(kiut) :: zdifcase ! case of the equation resolution (1->4) 125 REAL(wp), DIMENSION(kiut) :: zftrice ! solar radiation transmitted through the ice 126 REAL(wp), DIMENSION(kiut) :: zihic, zhsu 127 REAL(wp), DIMENSION(kiut,0:nlay_i) :: ztcond_i ! Ice thermal conductivity 128 REAL(wp), DIMENSION(kiut,0:nlay_i) :: zradtr_i ! Radiation transmitted through the ice 129 REAL(wp), DIMENSION(kiut,0:nlay_i) :: zradab_i ! Radiation absorbed in the ice 130 REAL(wp), DIMENSION(kiut,0:nlay_i) :: zkappa_i ! Kappa factor in the ice 131 REAL(wp), DIMENSION(kiut,0:nlay_i) :: ztiold ! Old temperature in the ice 132 REAL(wp), DIMENSION(kiut,0:nlay_i) :: zeta_i ! Eta factor in the ice 133 REAL(wp), DIMENSION(kiut,0:nlay_i) :: ztitemp ! Temporary temperature in the ice to check the convergence 134 REAL(wp), DIMENSION(kiut,0:nlay_i) :: zspeche_i ! Ice specific heat 135 REAL(wp), DIMENSION(kiut,0:nlay_i) :: z_i ! Vertical cotes of the layers in the ice 136 REAL(wp), DIMENSION(kiut,0:nlay_s) :: zradtr_s ! Radiation transmited through the snow 137 REAL(wp), DIMENSION(kiut,0:nlay_s) :: zradab_s ! Radiation absorbed in the snow 138 REAL(wp), DIMENSION(kiut,0:nlay_s) :: zkappa_s ! Kappa factor in the snow 139 REAL(wp), DIMENSION(kiut,0:nlay_s) :: zeta_s ! Eta factor in the snow 140 REAL(wp), DIMENSION(kiut,0:nlay_s) :: ztstemp ! Temporary temperature in the snow to check the convergence 141 REAL(wp), DIMENSION(kiut,0:nlay_s) :: ztsold ! Temporary temperature in the snow 142 REAL(wp), DIMENSION(kiut,0:nlay_s) :: z_s ! Vertical cotes of the layers in the snow 143 REAL(wp), DIMENSION(kiut,jkmax+2) :: zindterm ! Independent term 144 REAL(wp), DIMENSION(kiut,jkmax+2) :: zindtbis ! temporary independent term 145 REAL(wp), DIMENSION(kiut,jkmax+2) :: zdiagbis 146 REAL(wp), DIMENSION(kiut,jkmax+2,3) :: ztrid ! tridiagonal system terms 115 REAL(wp), POINTER, DIMENSION(:) :: ztfs ! ice melting point 116 REAL(wp), POINTER, DIMENSION(:) :: ztsuold ! old surface temperature (before the iterative procedure ) 117 REAL(wp), POINTER, DIMENSION(:) :: ztsuoldit ! surface temperature at previous iteration 118 REAL(wp), POINTER, DIMENSION(:) :: zh_i ! ice layer thickness 119 REAL(wp), POINTER, DIMENSION(:) :: zh_s ! snow layer thickness 120 REAL(wp), POINTER, DIMENSION(:) :: zfsw ! solar radiation absorbed at the surface 121 REAL(wp), POINTER, DIMENSION(:) :: zf ! surface flux function 122 REAL(wp), POINTER, DIMENSION(:) :: dzf ! derivative of the surface flux function 123 REAL(wp), POINTER, DIMENSION(:) :: zerrit ! current error on temperature 124 REAL(wp), POINTER, DIMENSION(:) :: zdifcase ! case of the equation resolution (1->4) 125 REAL(wp), POINTER, DIMENSION(:) :: zftrice ! solar radiation transmitted through the ice 126 REAL(wp), POINTER, DIMENSION(:) :: zihic, zhsu 127 REAL(wp), POINTER, DIMENSION(:,:) :: ztcond_i ! Ice thermal conductivity 128 REAL(wp), POINTER, DIMENSION(:,:) :: zradtr_i ! Radiation transmitted through the ice 129 REAL(wp), POINTER, DIMENSION(:,:) :: zradab_i ! Radiation absorbed in the ice 130 REAL(wp), POINTER, DIMENSION(:,:) :: zkappa_i ! Kappa factor in the ice 131 REAL(wp), POINTER, DIMENSION(:,:) :: ztiold ! Old temperature in the ice 132 REAL(wp), POINTER, DIMENSION(:,:) :: zeta_i ! Eta factor in the ice 133 REAL(wp), POINTER, DIMENSION(:,:) :: ztitemp ! Temporary temperature in the ice to check the convergence 134 REAL(wp), POINTER, DIMENSION(:,:) :: zspeche_i ! Ice specific heat 135 REAL(wp), POINTER, DIMENSION(:,:) :: z_i ! Vertical cotes of the layers in the ice 136 REAL(wp), POINTER, DIMENSION(:,:) :: zradtr_s ! Radiation transmited through the snow 137 REAL(wp), POINTER, DIMENSION(:,:) :: zradab_s ! Radiation absorbed in the snow 138 REAL(wp), POINTER, DIMENSION(:,:) :: zkappa_s ! Kappa factor in the snow 139 REAL(wp), POINTER, DIMENSION(:,:) :: zeta_s ! Eta factor in the snow 140 REAL(wp), POINTER, DIMENSION(:,:) :: ztstemp ! Temporary temperature in the snow to check the convergence 141 REAL(wp), POINTER, DIMENSION(:,:) :: ztsold ! Temporary temperature in the snow 142 REAL(wp), POINTER, DIMENSION(:,:) :: z_s ! Vertical cotes of the layers in the snow 143 REAL(wp), POINTER, DIMENSION(:,:) :: zindterm ! Independent term 144 REAL(wp), POINTER, DIMENSION(:,:) :: zindtbis ! temporary independent term 145 REAL(wp), POINTER, DIMENSION(:,:) :: zdiagbis 146 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrid ! tridiagonal system terms 147 ! diag errors on heat 148 REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini 149 REAL(wp) :: zhfx_err 147 150 !!------------------------------------------------------------------ 148 151 ! 152 CALL wrk_alloc( jpij, numeqmin, numeqmax, isnow ) 153 CALL wrk_alloc( jpij, ztfs, ztsuold, ztsuoldit, zh_i, zh_s, zfsw ) 154 CALL wrk_alloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zhsu ) 155 CALL wrk_alloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztiold, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0) 156 CALL wrk_alloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsold, zeta_s, ztstemp, z_s, kjstart=0) 157 CALL wrk_alloc( jpij, jkmax+2, zindterm, zindtbis, zdiagbis ) 158 CALL wrk_alloc( jpij, jkmax+2, 3, ztrid ) 159 160 CALL wrk_alloc( jpij, zdq, zq_ini ) 161 162 ! --- diag error on heat diffusion - PART 1 --- ! 163 zdq(:) = 0._wp ; zq_ini(:) = 0._wp 164 DO ji = kideb, kiut 165 zq_ini(ji) = ( SUM( q_i_b(ji,1:nlay_i) ) * ht_i_b(ji) / REAL( nlay_i ) + & 166 & SUM( q_s_b(ji,1:nlay_s) ) * ht_s_b(ji) / REAL( nlay_s ) ) 167 END DO 168 149 169 !------------------------------------------------------------------------------! 150 170 ! 1) Initialization ! 151 171 !------------------------------------------------------------------------------! 152 ! 172 ! clem clean: replace just ztfs by rtt 153 173 DO ji = kideb , kiut 154 174 ! is there snow or not 155 175 isnow(ji)= NINT( 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_b(ji) ) ) ) 156 176 ! surface temperature of fusion 157 !!gm ??? ztfs(ji) = rtt !!!????158 177 ztfs(ji) = REAL( isnow(ji) ) * rtt + REAL( 1 - isnow(ji) ) * rtt 159 178 ! layer thickness … … 194 213 ! zfsw = (1-i0).qsr_ice is absorbed at the surface 195 214 ! zftrice = io.qsr_ice is below the surface 196 ! f stbif= io.qsr_ice.exp(-k(h_i)) transmitted below the ice215 ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice 197 216 198 217 DO ji = kideb , kiut … … 253 272 254 273 DO ji = kideb, kiut ! Radiation transmitted below the ice 255 fstbif_1d(ji) = fstbif_1d(ji) + iatte_1d(ji) * zradtr_i(ji,nlay_i) * a_i_b(ji) / at_i_b(ji) ! clem modif 256 END DO 257 258 ! +++++ 259 ! just to check energy conservation 260 DO ji = kideb, kiut 261 ii = MOD( npb(ji) - 1 , jpi ) + 1 262 ij = ( npb(ji) - 1 ) / jpi + 1 263 fstroc(ii,ij,jl) = iatte_1d(ji) * zradtr_i(ji,nlay_i) ! clem modif 264 END DO 265 ! +++++ 266 267 DO layer = 1, nlay_i 268 DO ji = kideb, kiut 269 radab(ji,layer) = zradab_i(ji,layer) 270 END DO 274 !!!ftr_ice_1d(ji) = ftr_ice_1d(ji) + iatte_1d(ji) * zradtr_i(ji,nlay_i) * a_i_b(ji) / at_i_b(ji) ! clem modif 275 ftr_ice_1d(ji) = zradtr_i(ji,nlay_i) 271 276 END DO 272 277 … … 279 284 ztsuold (ji) = t_su_b(ji) ! temperature at the beg of iter pr. 280 285 ztsuoldit(ji) = t_su_b(ji) ! temperature at the previous iter 281 t_su_b (ji) = MIN( t_su_b(ji), ztfs(ji) -0.00001 )! necessary286 t_su_b (ji) = MIN( t_su_b(ji), ztfs(ji) - ztsu_err ) ! necessary 282 287 zerrit (ji) = 1000._wp ! initial value of error 283 288 END DO … … 403 408 ! 404 409 DO ji = kideb , kiut 405 406 410 ! update of the non solar flux according to the update in T_su 407 qnsr_ice_1d(ji) = qnsr_ice_1d(ji) + dqns_ice_1d(ji) * & 408 ( t_su_b(ji) - ztsuoldit(ji) ) 411 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_b(ji) - ztsuoldit(ji) ) 409 412 410 413 ! update incoming flux 411 414 zf(ji) = zfsw(ji) & ! net absorbed solar radiation 412 + qns r_ice_1d(ji)! non solar total flux415 + qns_ice_1d(ji) ! non solar total flux 413 416 ! (LWup, LWdw, SH, LH) 414 415 417 END DO 416 418 … … 678 680 DO layer = 1, nlay_s 679 681 DO ji = kideb , kiut 680 ii = MOD( npb(ji) - 1, jpi ) + 1681 ij = ( npb(ji) - 1 ) / jpi + 1682 682 t_s_b(ji,layer) = MAX( MIN( t_s_b(ji,layer), rtt ), 190._wp ) 683 683 zerrit(ji) = MAX(zerrit(ji),ABS(t_s_b(ji,layer) - ztstemp(ji,layer))) … … 713 713 !-------------------------------------------------------------------------! 714 714 DO ji = kideb, kiut 715 #if ! defined key_coupled716 715 ! forced mode only : update of latent heat fluxes (sublimation) (always >=0, upward flux) 717 qla_ice_1d (ji) = MAX( 0._wp, qla_ice_1d (ji) + dqla_ice_1d(ji) * ( t_su_b(ji) - ztsuold(ji) ) ) 718 #endif 716 IF( .NOT. lk_cpl) qla_ice_1d (ji) = MAX( 0._wp, qla_ice_1d (ji) + dqla_ice_1d(ji) * ( t_su_b(ji) - ztsuold(ji) ) ) 719 717 ! ! surface ice conduction flux 720 718 isnow(ji) = NINT( 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_b(ji) ) ) ) … … 725 723 END DO 726 724 727 !-------------------------! 728 ! Heat conservation ! 729 !-------------------------! 730 IF( con_i .AND. jiindex_1d > 0 ) THEN 725 !----------------------------------------- 726 ! Heat flux used to warm/cool ice in W.m-2 727 !----------------------------------------- 728 DO ji = kideb, kiut 729 IF( t_su_b(ji) < rtt ) THEN ! case T_su < 0degC 730 hfx_dif_1d(ji) = hfx_dif_1d(ji) + ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_b(ji) 731 ELSE ! case T_su = 0degC 732 hfx_dif_1d(ji) = hfx_dif_1d(ji) + ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_b(ji) 733 ENDIF 734 END DO 735 736 ! --- computes sea ice energy of melting compulsory for limthd_dh --- ! 737 CALL lim_thd_enmelt( kideb, kiut ) 738 739 ! --- diag error on heat diffusion - PART 2 --- ! 740 DO ji = kideb, kiut 741 zdq(ji) = - zq_ini(ji) + ( SUM( q_i_b(ji,1:nlay_i) ) * ht_i_b(ji) / REAL( nlay_i ) + & 742 & SUM( q_s_b(ji,1:nlay_s) ) * ht_s_b(ji) / REAL( nlay_s ) ) 743 zhfx_err = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice ) 744 hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err * a_i_b(ji) 745 ! --- correction of qns_ice and surface conduction flux --- ! 746 qns_ice_1d(ji) = qns_ice_1d(ji) - zhfx_err 747 fc_su (ji) = fc_su (ji) - zhfx_err 748 ! --- Heat flux at the ice surface in W.m-2 --- ! 749 ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 750 hfx_in (ii,ij) = hfx_in (ii,ij) + a_i_b(ji) * ( qsr_ice_1d(ji) + qns_ice_1d(ji) ) 751 END DO 752 753 ! 754 CALL wrk_dealloc( jpij, numeqmin, numeqmax, isnow ) 755 CALL wrk_dealloc( jpij, ztfs, ztsuold, ztsuoldit, zh_i, zh_s, zfsw ) 756 CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zhsu ) 757 CALL wrk_dealloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztiold, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 ) 758 CALL wrk_dealloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsold, zeta_s, ztstemp, z_s, kjstart = 0 ) 759 CALL wrk_dealloc( jpij, jkmax+2, zindterm, zindtbis, zdiagbis ) 760 CALL wrk_dealloc( jpij, jkmax+2, 3, ztrid ) 761 CALL wrk_dealloc( jpij, zdq, zq_ini ) 762 763 END SUBROUTINE lim_thd_dif 764 765 SUBROUTINE lim_thd_enmelt( kideb, kiut ) 766 !!----------------------------------------------------------------------- 767 !! *** ROUTINE lim_thd_enmelt *** 768 !! 769 !! ** Purpose : Computes sea ice energy of melting q_i (J.m-3) from temperature 770 !! 771 !! ** Method : Formula (Bitz and Lipscomb, 1999) 772 !!------------------------------------------------------------------- 773 INTEGER, INTENT(in) :: kideb, kiut ! bounds for the spatial loop 774 ! 775 INTEGER :: ji, jk ! dummy loop indices 776 REAL(wp) :: ztmelts, zindb ! local scalar 777 !!------------------------------------------------------------------- 778 ! 779 DO jk = 1, nlay_i ! Sea ice energy of melting 731 780 DO ji = kideb, kiut 732 ! Upper snow value 733 fc_s(ji,0) = - REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * ( t_s_b(ji,1) - t_su_b(ji) ) 734 ! Bott. snow value 735 fc_s(ji,1) = - REAL( isnow(ji) ) * zkappa_s(ji,1) * ( t_i_b(ji,1) - t_s_b(ji,1) ) 736 END DO 737 DO ji = kideb, kiut ! Upper ice layer 738 fc_i(ji,0) = - REAL( isnow(ji) ) * & ! interface flux if there is snow 739 ( zkappa_i(ji,0) * ( t_i_b(ji,1) - t_s_b(ji,nlay_s ) ) ) & 740 - REAL( 1 - isnow(ji) ) * ( zkappa_i(ji,0) * & 741 zg1 * ( t_i_b(ji,1) - t_su_b(ji) ) ) ! upper flux if not 742 END DO 743 DO layer = 1, nlay_i - 1 ! Internal ice layers 744 DO ji = kideb, kiut 745 fc_i(ji,layer) = - zkappa_i(ji,layer) * ( t_i_b(ji,layer+1) - t_i_b(ji,layer) ) 746 ii = MOD( npb(ji) - 1, jpi ) + 1 747 ij = ( npb(ji) - 1 ) / jpi + 1 748 END DO 749 END DO 750 DO ji = kideb, kiut ! Bottom ice layers 751 fc_i(ji,nlay_i) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) ) 752 END DO 753 ENDIF 781 ztmelts = - tmut * s_i_b(ji,jk) + rtt 782 zindb = MAX( 0._wp , SIGN( 1._wp , -(t_i_b(ji,jk) - rtt) - epsi10 ) ) 783 q_i_b(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_b(ji,jk) ) & 784 & + lfus * ( 1.0 - zindb * ( ztmelts-rtt ) / MIN( t_i_b(ji,jk)-rtt, -epsi10 ) ) & 785 & - rcp * ( ztmelts-rtt ) ) 786 END DO 787 END DO 788 DO jk = 1, nlay_s ! Snow energy of melting 789 DO ji = kideb, kiut 790 q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) 791 END DO 792 END DO 754 793 ! 755 END SUBROUTINE lim_thd_ dif794 END SUBROUTINE lim_thd_enmelt 756 795 757 796 #else
Note: See TracChangeset
for help on using the changeset viewer.
