- Timestamp:
- 2014-05-12T22:46:18+02:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r4332 r4634 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 … … 111 112 REAL(wp) :: zraext_s = 1.e+8_wp ! extinction coefficient of radiation in the snow 112 113 REAL(wp) :: zkimin = 0.10_wp ! minimum ice thermal conductivity 114 REAL(wp) :: ztsu_err = 1.e-5_wp ! range around which t_su is considered as 0°C 113 115 REAL(wp) :: ztmelt_i ! ice melting temperature 114 116 REAL(wp) :: zerritmax ! current maximal error on temperature … … 145 147 REAL(wp), DIMENSION(kiut,jkmax+2) :: zdiagbis 146 148 REAL(wp), DIMENSION(kiut,jkmax+2,3) :: ztrid ! tridiagonal system terms 149 REAL(wp) :: ztemp ! local scalar 147 150 !!------------------------------------------------------------------ 148 151 ! … … 150 153 ! 1) Initialization ! 151 154 !------------------------------------------------------------------------------! 152 ! 155 ! clem clean: replace just ztfs by rtt 153 156 DO ji = kideb , kiut 154 157 ! is there snow or not 155 158 isnow(ji)= NINT( 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_b(ji) ) ) ) 156 159 ! surface temperature of fusion 157 !!gm ??? ztfs(ji) = rtt !!!????158 160 ztfs(ji) = REAL( isnow(ji) ) * rtt + REAL( 1 - isnow(ji) ) * rtt 159 161 ! layer thickness … … 194 196 ! zfsw = (1-i0).qsr_ice is absorbed at the surface 195 197 ! zftrice = io.qsr_ice is below the surface 196 ! f stbif= io.qsr_ice.exp(-k(h_i)) transmitted below the ice198 ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice 197 199 198 200 DO ji = kideb , kiut … … 253 255 254 256 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 257 !!!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 258 ftr_ice_1d(ji) = zradtr_i(ji,nlay_i) 271 259 END DO 272 260 … … 279 267 ztsuold (ji) = t_su_b(ji) ! temperature at the beg of iter pr. 280 268 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 ) ! necessary 269 t_su_b (ji) = MIN( t_su_b(ji), ztfs(ji) - ztsu_err ) ! necessary 270 !!ztsuold (ji) = t_su_b(ji) ! temperature at the beg of iter pr. 271 !!ztsuoldit(ji) = t_su_b(ji) ! temperature at the previous iter 282 272 zerrit (ji) = 1000._wp ! initial value of error 283 273 END DO … … 328 318 DO layer = 1, nlay_i-1 329 319 DO ji = kideb , kiut 330 ztcond_i(ji,layer) = rcdic + 0.090_wp * ( s_i_b(ji,layer) + s_i_b(ji,layer+1) ) & 331 & / MIN(-2.0_wp * epsi10, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt) & 332 & - 0.0055_wp* ( t_i_b(ji,layer) + t_i_b(ji,layer+1) - 2.0*rtt ) 320 ztemp = t_i_b(ji,layer) + t_i_b(ji,layer+1) - 2._wp * rtt 321 ztcond_i(ji,layer) = rcdic + 0.0900_wp * ( s_i_b(ji,layer) + s_i_b(ji,layer+1) ) & 322 & / MIN( -2.0_wp * epsi10, ztemp ) & 323 & - 0.0055_wp * ztemp 333 324 ztcond_i(ji,layer) = MAX( ztcond_i(ji,layer), zkimin ) 334 325 END DO 335 326 END DO 336 327 DO ji = kideb , kiut 337 ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_b(ji,nlay_i) / MIN(-epsi10,t_bo_b(ji)-rtt) & 338 & - 0.011_wp * ( t_bo_b(ji) - rtt ) 328 ztemp = t_bo_b(ji) - rtt 329 ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_b(ji,nlay_i) / MIN( -epsi10, ztemp ) & 330 & - 0.011_wp * ztemp 339 331 ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin ) 340 332 END DO … … 405 397 406 398 ! 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) ) 399 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_b(ji) - ztsuoldit(ji) ) 409 400 410 401 ! update incoming flux 411 402 zf(ji) = zfsw(ji) & ! net absorbed solar radiation 412 + qns r_ice_1d(ji)! non solar total flux403 + qns_ice_1d(ji) ! non solar total flux 413 404 ! (LWup, LWdw, SH, LH) 414 405 406 ! heat flux used to change surface temperature 407 !hfx_tot_1d(ji) = hfx_tot_1d(ji) + dqns_ice_1d(ji) * ( t_su_b(ji) - ztsuoldit(ji) ) * a_i_b(ji) 415 408 END DO 416 409 … … 713 706 !-------------------------------------------------------------------------! 714 707 DO ji = kideb, kiut 715 #if ! defined key_coupled716 708 ! 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 709 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 710 ! ! surface ice conduction flux 720 711 isnow(ji) = NINT( 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_b(ji) ) ) ) … … 725 716 END DO 726 717 727 !-------------------------! 728 ! Heat conservation ! 729 !-------------------------! 730 IF( con_i .AND. jiindex_1d > 0 ) THEN 731 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 718 !----------------------------------------- 719 ! Heat flux used to warm/cool ice in W.m-2 720 !----------------------------------------- 721 DO ji = kideb, kiut 722 IF( t_su_b(ji) < rtt ) THEN ! case T_su < 0degC 723 hfx_tot_1d(ji) = hfx_tot_1d(ji) + ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_b(ji) 724 ELSE ! case T_su = 0degC 725 hfx_tot_1d(ji) = hfx_tot_1d(ji) + ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_b(ji) 726 ENDIF 727 END DO 728 754 729 ! 755 730 END SUBROUTINE lim_thd_dif
Note: See TracChangeset
for help on using the changeset viewer.
