- Timestamp:
- 2015-02-11T11:50:34+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO7_STARTHOUR/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r4333 r5075 25 25 USE wrk_nemo ! work arrays 26 26 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 27 USE sbc_oce, 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 34 !!---------------------------------------------------------------------- 35 35 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) … … 39 39 CONTAINS 40 40 41 SUBROUTINE lim_thd_dif( kideb , kiut , jl)41 SUBROUTINE lim_thd_dif( kideb , kiut ) 42 42 !!------------------------------------------------------------------ 43 43 !! *** ROUTINE lim_thd_dif *** … … 74 74 !! 75 75 !! ** Inputs / Ouputs : (global commons) 76 !! surface temperature : t_su_ b77 !! ice/snow temperatures : t_i_ b, t_s_b78 !! ice salinities : s_i_ b76 !! surface temperature : t_su_1d 77 !! ice/snow temperatures : t_i_1d, t_s_1d 78 !! ice salinities : s_i_1d 79 79 !! number of layers in the ice/snow: nlay_i, nlay_s 80 80 !! profile of the ice/snow layers : z_i, z_s 81 !! total ice/snow thickness : ht_i_ b, ht_s_b81 !! total ice/snow thickness : ht_i_1d, ht_s_1d 82 82 !! 83 83 !! ** External : … … 91 91 !! (04-2007) Energy conservation tested by M. Vancoppenolle 92 92 !!------------------------------------------------------------------ 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 93 INTEGER , INTENT(in) :: kideb, kiut ! Start/End point on which the the computation is applied 96 94 97 95 !! * Local variables … … 99 97 INTEGER :: ii, ij ! temporary dummy loop index 100 98 INTEGER :: numeq ! current reference number of equation 101 INTEGER :: layer! vertical dummy loop index99 INTEGER :: jk ! vertical dummy loop index 102 100 INTEGER :: nconv ! number of iterations in iterative procedure 103 101 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 snow102 INTEGER, POINTER, DIMENSION(:) :: numeqmin ! reference number of top equation 103 INTEGER, POINTER, DIMENSION(:) :: numeqmax ! reference number of bottom equation 104 INTEGER, POINTER, DIMENSION(:) :: isnow ! switch for presence (1) or absence (0) of snow 107 105 REAL(wp) :: zg1s = 2._wp ! for the tridiagonal system 108 106 REAL(wp) :: zg1 = 2._wp ! 109 107 REAL(wp) :: zgamma = 18009._wp ! for specific heat 110 108 REAL(wp) :: zbeta = 0.117_wp ! for thermal conductivity (could be 0.13) 111 REAL(wp) :: zraext_s = 1 .e+8_wp! extinction coefficient of radiation in the snow109 REAL(wp) :: zraext_s = 10._wp ! extinction coefficient of radiation in the snow 112 110 REAL(wp) :: zkimin = 0.10_wp ! minimum ice thermal conductivity 111 REAL(wp) :: ztsu_err = 1.e-5_wp ! range around which t_su is considered as 0°C 113 112 REAL(wp) :: ztmelt_i ! ice melting temperature 114 113 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 114 REAL(wp), POINTER, DIMENSION(:) :: ztfs ! ice melting point 115 REAL(wp), POINTER, DIMENSION(:) :: ztsub ! old surface temperature (before the iterative procedure ) 116 REAL(wp), POINTER, DIMENSION(:) :: ztsubit ! surface temperature at previous iteration 117 REAL(wp), POINTER, DIMENSION(:) :: zh_i ! ice layer thickness 118 REAL(wp), POINTER, DIMENSION(:) :: zh_s ! snow layer thickness 119 REAL(wp), POINTER, DIMENSION(:) :: zfsw ! solar radiation absorbed at the surface 120 REAL(wp), POINTER, DIMENSION(:) :: zf ! surface flux function 121 REAL(wp), POINTER, DIMENSION(:) :: dzf ! derivative of the surface flux function 122 REAL(wp), POINTER, DIMENSION(:) :: zerrit ! current error on temperature 123 REAL(wp), POINTER, DIMENSION(:) :: zdifcase ! case of the equation resolution (1->4) 124 REAL(wp), POINTER, DIMENSION(:) :: zftrice ! solar radiation transmitted through the ice 125 REAL(wp), POINTER, DIMENSION(:) :: zihic, zhsu 126 REAL(wp), POINTER, DIMENSION(:,:) :: ztcond_i ! Ice thermal conductivity 127 REAL(wp), POINTER, DIMENSION(:,:) :: zradtr_i ! Radiation transmitted through the ice 128 REAL(wp), POINTER, DIMENSION(:,:) :: zradab_i ! Radiation absorbed in the ice 129 REAL(wp), POINTER, DIMENSION(:,:) :: zkappa_i ! Kappa factor in the ice 130 REAL(wp), POINTER, DIMENSION(:,:) :: ztib ! Old temperature in the ice 131 REAL(wp), POINTER, DIMENSION(:,:) :: zeta_i ! Eta factor in the ice 132 REAL(wp), POINTER, DIMENSION(:,:) :: ztitemp ! Temporary temperature in the ice to check the convergence 133 REAL(wp), POINTER, DIMENSION(:,:) :: zspeche_i ! Ice specific heat 134 REAL(wp), POINTER, DIMENSION(:,:) :: z_i ! Vertical cotes of the layers in the ice 135 REAL(wp), POINTER, DIMENSION(:,:) :: zradtr_s ! Radiation transmited through the snow 136 REAL(wp), POINTER, DIMENSION(:,:) :: zradab_s ! Radiation absorbed in the snow 137 REAL(wp), POINTER, DIMENSION(:,:) :: zkappa_s ! Kappa factor in the snow 138 REAL(wp), POINTER, DIMENSION(:,:) :: zeta_s ! Eta factor in the snow 139 REAL(wp), POINTER, DIMENSION(:,:) :: ztstemp ! Temporary temperature in the snow to check the convergence 140 REAL(wp), POINTER, DIMENSION(:,:) :: ztsb ! Temporary temperature in the snow 141 REAL(wp), POINTER, DIMENSION(:,:) :: z_s ! Vertical cotes of the layers in the snow 142 REAL(wp), POINTER, DIMENSION(:,:) :: zswiterm ! Independent term 143 REAL(wp), POINTER, DIMENSION(:,:) :: zswitbis ! temporary independent term 144 REAL(wp), POINTER, DIMENSION(:,:) :: zdiagbis 145 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrid ! tridiagonal system terms 146 ! diag errors on heat 147 REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini, zhfx_err 147 148 !!------------------------------------------------------------------ 148 149 ! 150 CALL wrk_alloc( jpij, numeqmin, numeqmax, isnow ) 151 CALL wrk_alloc( jpij, ztfs, ztsub, ztsubit, zh_i, zh_s, zfsw ) 152 CALL wrk_alloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zhsu ) 153 CALL wrk_alloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0) 154 CALL wrk_alloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart=0) 155 CALL wrk_alloc( jpij, nlay_i+3, zswiterm, zswitbis, zdiagbis ) 156 CALL wrk_alloc( jpij, nlay_i+3, 3, ztrid ) 157 158 CALL wrk_alloc( jpij, zdq, zq_ini, zhfx_err ) 159 160 ! --- diag error on heat diffusion - PART 1 --- ! 161 zdq(:) = 0._wp ; zq_ini(:) = 0._wp 162 DO ji = kideb, kiut 163 zq_ini(ji) = ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) / REAL( nlay_i ) + & 164 & SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) / REAL( nlay_s ) ) 165 END DO 166 149 167 !------------------------------------------------------------------------------! 150 168 ! 1) Initialization ! 151 169 !------------------------------------------------------------------------------! 152 ! 170 ! clem clean: replace just ztfs by rtt 153 171 DO ji = kideb , kiut 154 172 ! is there snow or not 155 isnow(ji)= NINT( 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_ b(ji) ) ) )173 isnow(ji)= NINT( 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) ) ) 156 174 ! surface temperature of fusion 157 !!gm ??? ztfs(ji) = rtt !!!????158 175 ztfs(ji) = REAL( isnow(ji) ) * rtt + REAL( 1 - isnow(ji) ) * rtt 159 176 ! layer thickness 160 zh_i(ji) = ht_i_ b(ji) / REAL( nlay_i )161 zh_s(ji) = ht_s_ b(ji) / REAL( nlay_s )177 zh_i(ji) = ht_i_1d(ji) / REAL( nlay_i ) 178 zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s ) 162 179 END DO 163 180 … … 169 186 z_i(:,0) = 0._wp ! vert. coord. of the up. lim. of the 1st ice layer 170 187 171 DO layer= 1, nlay_s ! vert. coord of the up. lim. of the layer-th snow layer172 DO ji = kideb , kiut 173 z_s(ji, layer) = z_s(ji,layer-1) + ht_s_b(ji) / REAL( nlay_s )174 END DO 175 END DO 176 177 DO layer= 1, nlay_i ! vert. coord of the up. lim. of the layer-th ice layer178 DO ji = kideb , kiut 179 z_i(ji, layer) = z_i(ji,layer-1) + ht_i_b(ji) / REAL( nlay_i )188 DO jk = 1, nlay_s ! vert. coord of the up. lim. of the layer-th snow layer 189 DO ji = kideb , kiut 190 z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) / REAL( nlay_s ) 191 END DO 192 END DO 193 194 DO jk = 1, nlay_i ! vert. coord of the up. lim. of the layer-th ice layer 195 DO ji = kideb , kiut 196 z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) / REAL( nlay_i ) 180 197 END DO 181 198 END DO … … 194 211 ! zfsw = (1-i0).qsr_ice is absorbed at the surface 195 212 ! zftrice = io.qsr_ice is below the surface 196 ! f stbif= io.qsr_ice.exp(-k(h_i)) transmitted below the ice213 ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice 197 214 198 215 DO ji = kideb , kiut 199 216 ! switches 200 isnow(ji) = NINT( 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_ b(ji) ) ) )217 isnow(ji) = NINT( 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) ) ) 201 218 ! hs > 0, isnow = 1 202 219 zhsu (ji) = hnzst ! threshold for the computation of i0 203 zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_ b(ji) / zhsu(ji) ) )220 zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) / zhsu(ji) ) ) 204 221 205 222 i0(ji) = REAL( 1 - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) ) … … 208 225 ! a function of the cloud cover 209 226 ! 210 !i0(ji) = (1.0-FLOAT(isnow(ji)))*3.0/(100*ht_s_ b(ji)+10.0)227 !i0(ji) = (1.0-FLOAT(isnow(ji)))*3.0/(100*ht_s_1d(ji)+10.0) 211 228 !formula used in Cice 212 229 END DO … … 230 247 END DO 231 248 232 DO layer= 1, nlay_s ! Radiation through snow249 DO jk = 1, nlay_s ! Radiation through snow 233 250 DO ji = kideb, kiut 234 251 ! ! radiation transmitted below the layer-th snow layer 235 zradtr_s(ji, layer) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,layer) ) ) )252 zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,jk) ) ) ) 236 253 ! ! radiation absorbed by the layer-th snow layer 237 zradab_s(ji, layer) = zradtr_s(ji,layer-1) - zradtr_s(ji,layer)254 zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk) 238 255 END DO 239 256 END DO … … 243 260 END DO 244 261 245 DO layer= 1, nlay_i ! Radiation through ice262 DO jk = 1, nlay_i ! Radiation through ice 246 263 DO ji = kideb, kiut 247 264 ! ! radiation transmitted below the layer-th ice layer 248 zradtr_i(ji, layer) = zradtr_i(ji,0) * EXP( - kappa_i * ( MAX ( 0._wp , z_i(ji,layer) ) ) )265 zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - kappa_i * ( MAX ( 0._wp , z_i(ji,jk) ) ) ) 249 266 ! ! radiation absorbed by the layer-th ice layer 250 zradab_i(ji, layer) = zradtr_i(ji,layer-1) - zradtr_i(ji,layer)267 zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk) 251 268 END DO 252 269 END DO 253 270 254 271 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 272 ftr_ice_1d(ji) = zradtr_i(ji,nlay_i) 271 273 END DO 272 274 … … 277 279 ! 278 280 DO ji = kideb, kiut ! Old surface temperature 279 ztsu old (ji) = t_su_b(ji) ! temperature at the beg of iter pr.280 ztsu oldit(ji) = t_su_b(ji) ! temperature at the previous iter281 t_su_ b (ji) = MIN( t_su_b(ji), ztfs(ji)-0.00001 )! necessary281 ztsub (ji) = t_su_1d(ji) ! temperature at the beg of iter pr. 282 ztsubit(ji) = t_su_1d(ji) ! temperature at the previous iter 283 t_su_1d (ji) = MIN( t_su_1d(ji), ztfs(ji) - ztsu_err ) ! necessary 282 284 zerrit (ji) = 1000._wp ! initial value of error 283 285 END DO 284 286 285 DO layer= 1, nlay_s ! Old snow temperature286 DO ji = kideb , kiut 287 zts old(ji,layer) = t_s_b(ji,layer)288 END DO 289 END DO 290 291 DO layer= 1, nlay_i ! Old ice temperature292 DO ji = kideb , kiut 293 zti old(ji,layer) = t_i_b(ji,layer)287 DO jk = 1, nlay_s ! Old snow temperature 288 DO ji = kideb , kiut 289 ztsb(ji,jk) = t_s_1d(ji,jk) 290 END DO 291 END DO 292 293 DO jk = 1, nlay_i ! Old ice temperature 294 DO ji = kideb , kiut 295 ztib(ji,jk) = t_i_1d(ji,jk) 294 296 END DO 295 297 END DO … … 308 310 IF( thcon_i_swi == 0 ) THEN ! Untersteiner (1964) formula 309 311 DO ji = kideb , kiut 310 ztcond_i(ji,0) = rcdic + zbeta*s_i_ b(ji,1) / MIN(-epsi10,t_i_b(ji,1)-rtt)312 ztcond_i(ji,0) = rcdic + zbeta*s_i_1d(ji,1) / MIN(-epsi10,t_i_1d(ji,1)-rtt) 311 313 ztcond_i(ji,0) = MAX(ztcond_i(ji,0),zkimin) 312 314 END DO 313 DO layer= 1, nlay_i-1315 DO jk = 1, nlay_i-1 314 316 DO ji = kideb , kiut 315 ztcond_i(ji, layer) = rcdic + zbeta*( s_i_b(ji,layer) + s_i_b(ji,layer+1) ) / &316 MIN(-2.0_wp * epsi10, t_i_ b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt)317 ztcond_i(ji, layer) = MAX(ztcond_i(ji,layer),zkimin)317 ztcond_i(ji,jk) = rcdic + zbeta*( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) / & 318 MIN(-2.0_wp * epsi10, t_i_1d(ji,jk)+t_i_1d(ji,jk+1) - 2.0_wp * rtt) 319 ztcond_i(ji,jk) = MAX(ztcond_i(ji,jk),zkimin) 318 320 END DO 319 321 END DO … … 322 324 IF( thcon_i_swi == 1 ) THEN ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T 323 325 DO ji = kideb , kiut 324 ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_ b(ji,1) / MIN( -epsi10, t_i_b(ji,1)-rtt ) &325 & - 0.011_wp * ( t_i_ b(ji,1) - rtt )326 ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1)-rtt ) & 327 & - 0.011_wp * ( t_i_1d(ji,1) - rtt ) 326 328 ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin ) 327 329 END DO 328 DO layer= 1, nlay_i-1330 DO jk = 1, nlay_i-1 329 331 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 ) 333 ztcond_i(ji,layer) = MAX( ztcond_i(ji,layer), zkimin ) 332 ztcond_i(ji,jk) = rcdic + & 333 & 0.090_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) & 334 & / MIN(-2.0_wp * epsi10, t_i_1d(ji,jk)+t_i_1d(ji,jk+1) - 2.0_wp * rtt) & 335 & - 0.0055_wp* ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0*rtt ) 336 ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin ) 334 337 END DO 335 338 END DO 336 339 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 )340 ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN(-epsi10,t_bo_1d(ji)-rtt) & 341 & - 0.011_wp * ( t_bo_1d(ji) - rtt ) 339 342 ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin ) 340 343 END DO … … 352 355 END DO 353 356 354 DO layer= 1, nlay_s-1355 DO ji = kideb , kiut 356 zkappa_s(ji, layer) = 2.0 * rcdsn / &357 DO jk = 1, nlay_s-1 358 DO ji = kideb , kiut 359 zkappa_s(ji,jk) = 2.0 * rcdsn / & 357 360 MAX(epsi10,2.0*zh_s(ji)) 358 361 END DO 359 362 END DO 360 363 361 DO layer= 1, nlay_i-1364 DO jk = 1, nlay_i-1 362 365 DO ji = kideb , kiut 363 366 !-- Ice kappa factors 364 zkappa_i(ji, layer) = 2.0*ztcond_i(ji,layer)/ &367 zkappa_i(ji,jk) = 2.0*ztcond_i(ji,jk)/ & 365 368 MAX(epsi10,2.0*zh_i(ji)) 366 369 END DO … … 381 384 !------------------------------------------------------------------------------| 382 385 ! 383 DO layer= 1, nlay_i384 DO ji = kideb , kiut 385 ztitemp(ji, layer) = t_i_b(ji,layer)386 zspeche_i(ji, layer) = cpic + zgamma*s_i_b(ji,layer)/ &387 MAX((t_i_ b(ji,layer)-rtt)*(ztiold(ji,layer)-rtt),epsi10)388 zeta_i(ji, layer) = rdt_ice / MAX(rhoic*zspeche_i(ji,layer)*zh_i(ji), &386 DO jk = 1, nlay_i 387 DO ji = kideb , kiut 388 ztitemp(ji,jk) = t_i_1d(ji,jk) 389 zspeche_i(ji,jk) = cpic + zgamma*s_i_1d(ji,jk)/ & 390 MAX((t_i_1d(ji,jk)-rtt)*(ztib(ji,jk)-rtt),epsi10) 391 zeta_i(ji,jk) = rdt_ice / MAX(rhoic*zspeche_i(ji,jk)*zh_i(ji), & 389 392 epsi10) 390 393 END DO 391 394 END DO 392 395 393 DO layer= 1, nlay_s394 DO ji = kideb , kiut 395 ztstemp(ji, layer) = t_s_b(ji,layer)396 zeta_s(ji, layer) = rdt_ice / MAX(rhosn*cpic*zh_s(ji),epsi10)396 DO jk = 1, nlay_s 397 DO ji = kideb , kiut 398 ztstemp(ji,jk) = t_s_1d(ji,jk) 399 zeta_s(ji,jk) = rdt_ice / MAX(rhosn*cpic*zh_s(ji),epsi10) 397 400 END DO 398 401 END DO … … 402 405 !------------------------------------------------------------------------------| 403 406 ! 404 DO ji = kideb , kiut 405 406 ! 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) ) 409 407 IF( .NOT. lk_cpl ) THEN !--- forced atmosphere case 408 DO ji = kideb , kiut 409 ! update of the non solar flux according to the update in T_su 410 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsubit(ji) ) 411 END DO 412 ENDIF 413 414 ! Update incoming flux 415 DO ji = kideb , kiut 410 416 ! update incoming flux 411 417 zf(ji) = zfsw(ji) & ! net absorbed solar radiation 412 + qns r_ice_1d(ji)! non solar total flux418 + qns_ice_1d(ji) ! non solar total flux 413 419 ! (LWup, LWdw, SH, LH) 414 415 420 END DO 416 421 … … 427 432 !!ice interior terms (top equation has the same form as the others) 428 433 429 DO numeq=1, jkmax+2434 DO numeq=1,nlay_i+3 430 435 DO ji = kideb , kiut 431 436 ztrid(ji,numeq,1) = 0. 432 437 ztrid(ji,numeq,2) = 0. 433 438 ztrid(ji,numeq,3) = 0. 434 z indterm(ji,numeq)= 0.435 z indtbis(ji,numeq)= 0.439 zswiterm(ji,numeq)= 0. 440 zswitbis(ji,numeq)= 0. 436 441 zdiagbis(ji,numeq)= 0. 437 442 ENDDO … … 440 445 DO numeq = nlay_s + 2, nlay_s + nlay_i 441 446 DO ji = kideb , kiut 442 layer= numeq - nlay_s - 1443 ztrid(ji,numeq,1) = - zeta_i(ji, layer)*zkappa_i(ji,layer-1)444 ztrid(ji,numeq,2) = 1.0 + zeta_i(ji, layer)*(zkappa_i(ji,layer-1) + &445 zkappa_i(ji, layer))446 ztrid(ji,numeq,3) = - zeta_i(ji, layer)*zkappa_i(ji,layer)447 z indterm(ji,numeq) = ztiold(ji,layer) + zeta_i(ji,layer)* &448 zradab_i(ji, layer)447 jk = numeq - nlay_s - 1 448 ztrid(ji,numeq,1) = - zeta_i(ji,jk)*zkappa_i(ji,jk-1) 449 ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,jk)*(zkappa_i(ji,jk-1) + & 450 zkappa_i(ji,jk)) 451 ztrid(ji,numeq,3) = - zeta_i(ji,jk)*zkappa_i(ji,jk) 452 zswiterm(ji,numeq) = ztib(ji,jk) + zeta_i(ji,jk)* & 453 zradab_i(ji,jk) 449 454 END DO 450 455 ENDDO … … 457 462 + zkappa_i(ji,nlay_i-1) ) 458 463 ztrid(ji,numeq,3) = 0.0 459 z indterm(ji,numeq) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i)* &464 zswiterm(ji,numeq) = ztib(ji,nlay_i) + zeta_i(ji,nlay_i)* & 460 465 ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i)*zg1 & 461 * t_bo_ b(ji) )466 * t_bo_1d(ji) ) 462 467 ENDDO 463 468 464 469 465 470 DO ji = kideb , kiut 466 IF ( ht_s_ b(ji).gt.0.0 ) THEN471 IF ( ht_s_1d(ji).gt.0.0 ) THEN 467 472 ! 468 473 !------------------------------------------------------------------------------| … … 472 477 !!snow interior terms (bottom equation has the same form as the others) 473 478 DO numeq = 3, nlay_s + 1 474 layer= numeq - 1475 ztrid(ji,numeq,1) = - zeta_s(ji, layer)*zkappa_s(ji,layer-1)476 ztrid(ji,numeq,2) = 1.0 + zeta_s(ji, layer)*( zkappa_s(ji,layer-1) + &477 zkappa_s(ji, layer) )478 ztrid(ji,numeq,3) = - zeta_s(ji, layer)*zkappa_s(ji,layer)479 z indterm(ji,numeq) = ztsold(ji,layer) + zeta_s(ji,layer)* &480 zradab_s(ji, layer)479 jk = numeq - 1 480 ztrid(ji,numeq,1) = - zeta_s(ji,jk)*zkappa_s(ji,jk-1) 481 ztrid(ji,numeq,2) = 1.0 + zeta_s(ji,jk)*( zkappa_s(ji,jk-1) + & 482 zkappa_s(ji,jk) ) 483 ztrid(ji,numeq,3) = - zeta_s(ji,jk)*zkappa_s(ji,jk) 484 zswiterm(ji,numeq) = ztsb(ji,jk) + zeta_s(ji,jk)* & 485 zradab_s(ji,jk) 481 486 END DO 482 487 … … 484 489 IF ( nlay_i.eq.1 ) THEN 485 490 ztrid(ji,nlay_s+2,3) = 0.0 486 z indterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1)* &487 t_bo_ b(ji)491 zswiterm(ji,nlay_s+2) = zswiterm(ji,nlay_s+2) + zkappa_i(ji,1)* & 492 t_bo_1d(ji) 488 493 ENDIF 489 494 490 IF ( t_su_ b(ji) .LT. rtt ) THEN495 IF ( t_su_1d(ji) .LT. rtt ) THEN 491 496 492 497 !------------------------------------------------------------------------------| … … 501 506 ztrid(ji,1,2) = dzf(ji) - zg1s*zkappa_s(ji,0) 502 507 ztrid(ji,1,3) = zg1s*zkappa_s(ji,0) 503 z indterm(ji,1) = dzf(ji)*t_su_b(ji) - zf(ji)508 zswiterm(ji,1) = dzf(ji)*t_su_1d(ji) - zf(ji) 504 509 505 510 !!first layer of snow equation … … 507 512 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1)*(zkappa_s(ji,1) + zkappa_s(ji,0)*zg1s) 508 513 ztrid(ji,2,3) = - zeta_s(ji,1)* zkappa_s(ji,1) 509 z indterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1)*zradab_s(ji,1)514 zswiterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1)*zradab_s(ji,1) 510 515 511 516 ELSE … … 524 529 zkappa_s(ji,0) * zg1s ) 525 530 ztrid(ji,2,3) = - zeta_s(ji,1)*zkappa_s(ji,1) 526 z indterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * &531 zswiterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) * & 527 532 ( zradab_s(ji,1) + & 528 zkappa_s(ji,0) * zg1s * t_su_ b(ji) )533 zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 529 534 ENDIF 530 535 ELSE … … 534 539 !------------------------------------------------------------------------------| 535 540 ! 536 IF (t_su_ b(ji) .LT. rtt) THEN541 IF (t_su_1d(ji) .LT. rtt) THEN 537 542 ! 538 543 !------------------------------------------------------------------------------| … … 548 553 ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0)*zg1 549 554 ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0)*zg1 550 z indterm(ji,numeqmin(ji)) = dzf(ji)*t_su_b(ji) - zf(ji)555 zswiterm(ji,numeqmin(ji)) = dzf(ji)*t_su_1d(ji) - zf(ji) 551 556 552 557 !!first layer of ice equation … … 555 560 + zkappa_i(ji,0) * zg1 ) 556 561 ztrid(ji,numeqmin(ji)+1,3) = - zeta_i(ji,1)*zkappa_i(ji,1) 557 z indterm(ji,numeqmin(ji)+1)= ztiold(ji,1) + zeta_i(ji,1)*zradab_i(ji,1)562 zswiterm(ji,numeqmin(ji)+1)= ztib(ji,1) + zeta_i(ji,1)*zradab_i(ji,1) 558 563 559 564 !!case of only one layer in the ice (surface & ice equations are altered) … … 568 573 ztrid(ji,numeqmin(ji)+1,3) = 0.0 569 574 570 z indterm(ji,numeqmin(ji)+1) = ztiold(ji,1) + zeta_i(ji,1)* &571 ( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_ b(ji) )575 zswiterm(ji,numeqmin(ji)+1) = ztib(ji,1) + zeta_i(ji,1)* & 576 ( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji) ) 572 577 ENDIF 573 578 … … 588 593 zg1) 589 594 ztrid(ji,numeqmin(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1) 590 z indterm(ji,numeqmin(ji)) = ztiold(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + &591 zkappa_i(ji,0) * zg1 * t_su_ b(ji) )595 zswiterm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + & 596 zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 592 597 593 598 !!case of only one layer in the ice (surface & ice equations are altered) … … 597 602 zkappa_i(ji,1)) 598 603 ztrid(ji,numeqmin(ji),3) = 0.0 599 z indterm(ji,numeqmin(ji)) = ztiold(ji,1) + zeta_i(ji,1)* &600 (zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_ b(ji)) &601 + t_su_ b(ji)*zeta_i(ji,1)*zkappa_i(ji,0)*2.0604 zswiterm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1)* & 605 (zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji)) & 606 + t_su_1d(ji)*zeta_i(ji,1)*zkappa_i(ji,0)*2.0 602 607 ENDIF 603 608 … … 618 623 619 624 maxnumeqmax = 0 620 minnumeqmin = jkmax+4621 622 DO ji = kideb , kiut 623 z indtbis(ji,numeqmin(ji)) = zindterm(ji,numeqmin(ji))625 minnumeqmin = nlay_i+5 626 627 DO ji = kideb , kiut 628 zswitbis(ji,numeqmin(ji)) = zswiterm(ji,numeqmin(ji)) 624 629 zdiagbis(ji,numeqmin(ji)) = ztrid(ji,numeqmin(ji),2) 625 630 minnumeqmin = MIN(numeqmin(ji),minnumeqmin) … … 627 632 END DO 628 633 629 DO layer= minnumeqmin+1, maxnumeqmax630 DO ji = kideb , kiut 631 numeq = min(max(numeqmin(ji)+1, layer),numeqmax(ji))634 DO jk = minnumeqmin+1, maxnumeqmax 635 DO ji = kideb , kiut 636 numeq = min(max(numeqmin(ji)+1,jk),numeqmax(ji)) 632 637 zdiagbis(ji,numeq) = ztrid(ji,numeq,2) - ztrid(ji,numeq,1)* & 633 638 ztrid(ji,numeq-1,3)/zdiagbis(ji,numeq-1) 634 z indtbis(ji,numeq) = zindterm(ji,numeq) - ztrid(ji,numeq,1)* &635 z indtbis(ji,numeq-1)/zdiagbis(ji,numeq-1)639 zswitbis(ji,numeq) = zswiterm(ji,numeq) - ztrid(ji,numeq,1)* & 640 zswitbis(ji,numeq-1)/zdiagbis(ji,numeq-1) 636 641 END DO 637 642 END DO … … 639 644 DO ji = kideb , kiut 640 645 ! ice temperatures 641 t_i_ b(ji,nlay_i) = zindtbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji))646 t_i_1d(ji,nlay_i) = zswitbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji)) 642 647 END DO 643 648 644 649 DO numeq = nlay_i + nlay_s + 1, nlay_s + 2, -1 645 650 DO ji = kideb , kiut 646 layer= numeq - nlay_s - 1647 t_i_ b(ji,layer) = (zindtbis(ji,numeq) - ztrid(ji,numeq,3)* &648 t_i_ b(ji,layer+1))/zdiagbis(ji,numeq)651 jk = numeq - nlay_s - 1 652 t_i_1d(ji,jk) = (zswitbis(ji,numeq) - ztrid(ji,numeq,3)* & 653 t_i_1d(ji,jk+1))/zdiagbis(ji,numeq) 649 654 END DO 650 655 END DO … … 652 657 DO ji = kideb , kiut 653 658 ! snow temperatures 654 IF (ht_s_ b(ji).GT.0._wp) &655 t_s_ b(ji,nlay_s) = (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) &656 * t_i_ b(ji,1))/zdiagbis(ji,nlay_s+1) &657 * MAX(0.0,SIGN(1.0,ht_s_ b(ji)))659 IF (ht_s_1d(ji).GT.0._wp) & 660 t_s_1d(ji,nlay_s) = (zswitbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) & 661 * t_i_1d(ji,1))/zdiagbis(ji,nlay_s+1) & 662 * MAX(0.0,SIGN(1.0,ht_s_1d(ji))) 658 663 659 664 ! surface temperature 660 isnow(ji) = NINT( 1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_ b(ji) ) ) )661 ztsu oldit(ji) = t_su_b(ji)662 IF( t_su_ b(ji) < ztfs(ji) ) &663 t_su_ b(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_b(ji,1) &664 & + REAL( 1 - isnow(ji) )*t_i_ b(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))665 isnow(ji) = NINT( 1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_1d(ji) ) ) ) 666 ztsubit(ji) = t_su_1d(ji) 667 IF( t_su_1d(ji) < ztfs(ji) ) & 668 t_su_1d(ji) = ( zswitbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_1d(ji,1) & 669 & + REAL( 1 - isnow(ji) )*t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji)) 665 670 END DO 666 671 ! … … 672 677 ! zerrit(ji) is a measure of error, it has to be under maxer_i_thd 673 678 DO ji = kideb , kiut 674 t_su_b(ji) = MAX( MIN( t_su_b(ji) , ztfs(ji) ) , 190._wp ) 675 zerrit(ji) = ABS( t_su_b(ji) - ztsuoldit(ji) ) 676 END DO 677 678 DO layer = 1, nlay_s 679 DO ji = kideb , kiut 680 ii = MOD( npb(ji) - 1, jpi ) + 1 681 ij = ( npb(ji) - 1 ) / jpi + 1 682 t_s_b(ji,layer) = MAX( MIN( t_s_b(ji,layer), rtt ), 190._wp ) 683 zerrit(ji) = MAX(zerrit(ji),ABS(t_s_b(ji,layer) - ztstemp(ji,layer))) 684 END DO 685 END DO 686 687 DO layer = 1, nlay_i 688 DO ji = kideb , kiut 689 ztmelt_i = -tmut * s_i_b(ji,layer) + rtt 690 t_i_b(ji,layer) = MAX(MIN(t_i_b(ji,layer),ztmelt_i), 190._wp) 691 zerrit(ji) = MAX(zerrit(ji),ABS(t_i_b(ji,layer) - ztitemp(ji,layer))) 679 t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , ztfs(ji) ) , 190._wp ) 680 zerrit(ji) = ABS( t_su_1d(ji) - ztsubit(ji) ) 681 END DO 682 683 DO jk = 1, nlay_s 684 DO ji = kideb , kiut 685 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rtt ), 190._wp ) 686 zerrit(ji) = MAX(zerrit(ji),ABS(t_s_1d(ji,jk) - ztstemp(ji,jk))) 687 END DO 688 END DO 689 690 DO jk = 1, nlay_i 691 DO ji = kideb , kiut 692 ztmelt_i = -tmut * s_i_1d(ji,jk) + rtt 693 t_i_1d(ji,jk) = MAX(MIN(t_i_1d(ji,jk),ztmelt_i), 190._wp) 694 zerrit(ji) = MAX(zerrit(ji),ABS(t_i_1d(ji,jk) - ztitemp(ji,jk))) 692 695 END DO 693 696 END DO … … 713 716 !-------------------------------------------------------------------------! 714 717 DO ji = kideb, kiut 715 #if ! defined key_coupled716 718 ! 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 719 IF( .NOT. lk_cpl) qla_ice_1d (ji) = MAX( 0._wp, qla_ice_1d (ji) + dqla_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) ) 719 720 ! ! surface ice conduction flux 720 isnow(ji) = NINT( 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_ b(ji) ) ) )721 fc_su(ji) = - REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * (t_s_ b(ji,1) - t_su_b(ji)) &722 & - REAL( 1 - isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_ b(ji,1) - t_su_b(ji))721 isnow(ji) = NINT( 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_1d(ji) ) ) ) 722 fc_su(ji) = - REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji)) & 723 & - REAL( 1 - isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_1d(ji,1) - t_su_1d(ji)) 723 724 ! ! bottom ice conduction flux 724 fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) ) 725 END DO 726 727 !-------------------------! 728 ! Heat conservation ! 729 !-------------------------! 730 IF( con_i .AND. jiindex_1d > 0 ) THEN 725 fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) ) 726 END DO 727 728 !----------------------------------------- 729 ! Heat flux used to warm/cool ice in W.m-2 730 !----------------------------------------- 731 DO ji = kideb, kiut 732 IF( t_su_1d(ji) < rtt ) THEN ! case T_su < 0degC 733 hfx_dif_1d(ji) = hfx_dif_1d(ji) + & 734 & ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji) 735 ELSE ! case T_su = 0degC 736 hfx_dif_1d(ji) = hfx_dif_1d(ji) + & 737 & ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji) 738 ENDIF 739 END DO 740 741 ! --- computes sea ice energy of melting compulsory for limthd_dh --- ! 742 CALL lim_thd_enmelt( kideb, kiut ) 743 744 ! --- diag conservation imbalance on heat diffusion - PART 2 --- ! 745 DO ji = kideb, kiut 746 zdq(ji) = - zq_ini(ji) + ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) / REAL( nlay_i ) + & 747 & SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) / REAL( nlay_s ) ) 748 zhfx_err(ji) = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice ) 749 hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err(ji) * a_i_1d(ji) 750 END DO 751 752 ! diagnose external surface (forced case) or bottom (forced case) from heat conservation 753 IF( .NOT. lk_cpl ) THEN ! --- forced case: qns_ice and fc_su are diagnosed 754 ! 731 755 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 756 qns_ice_1d(ji) = qns_ice_1d(ji) - zhfx_err(ji) 757 fc_su (ji) = fc_su(ji) - zhfx_err(ji) 758 END DO 759 ! 760 ELSE ! --- coupled case: ocean turbulent heat flux is diagnosed 761 ! 762 DO ji = kideb, kiut 763 fhtur_1d (ji) = fhtur_1d(ji) - zhfx_err(ji) 764 END DO 765 ! 753 766 ENDIF 767 768 ! --- compute diagnostic net heat flux at the surface of the snow-ice system (W.m2) 769 DO ji = kideb, kiut 770 ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 771 hfx_in (ii,ij) = hfx_in (ii,ij) + a_i_1d(ji) * ( qsr_ice_1d(ji) + qns_ice_1d(ji) ) 772 END DO 773 754 774 ! 775 CALL wrk_dealloc( jpij, numeqmin, numeqmax, isnow ) 776 CALL wrk_dealloc( jpij, ztfs, ztsub, ztsubit, zh_i, zh_s, zfsw ) 777 CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zhsu ) 778 CALL wrk_dealloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, & 779 & ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 ) 780 CALL wrk_dealloc( jpij, nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart = 0 ) 781 CALL wrk_dealloc( jpij, nlay_i+3, zswiterm, zswitbis, zdiagbis ) 782 CALL wrk_dealloc( jpij, nlay_i+3, 3, ztrid ) 783 CALL wrk_dealloc( jpij, zdq, zq_ini, zhfx_err ) 784 755 785 END SUBROUTINE lim_thd_dif 786 787 SUBROUTINE lim_thd_enmelt( kideb, kiut ) 788 !!----------------------------------------------------------------------- 789 !! *** ROUTINE lim_thd_enmelt *** 790 !! 791 !! ** Purpose : Computes sea ice energy of melting q_i (J.m-3) from temperature 792 !! 793 !! ** Method : Formula (Bitz and Lipscomb, 1999) 794 !!------------------------------------------------------------------- 795 INTEGER, INTENT(in) :: kideb, kiut ! bounds for the spatial loop 796 ! 797 INTEGER :: ji, jk ! dummy loop indices 798 REAL(wp) :: ztmelts ! local scalar 799 !!------------------------------------------------------------------- 800 ! 801 DO jk = 1, nlay_i ! Sea ice energy of melting 802 DO ji = kideb, kiut 803 ztmelts = - tmut * s_i_1d(ji,jk) + rtt 804 rswitch = MAX( 0._wp , SIGN( 1._wp , -(t_i_1d(ji,jk) - rtt) - epsi10 ) ) 805 q_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) ) & 806 & + lfus * ( 1.0 - rswitch * ( ztmelts-rtt ) / MIN( t_i_1d(ji,jk)-rtt, -epsi10 ) ) & 807 & - rcp * ( ztmelts-rtt ) ) 808 END DO 809 END DO 810 DO jk = 1, nlay_s ! Snow energy of melting 811 DO ji = kideb, kiut 812 q_s_1d(ji,jk) = rhosn * ( cpic * ( rtt - t_s_1d(ji,jk) ) + lfus ) 813 END DO 814 END DO 815 ! 816 END SUBROUTINE lim_thd_enmelt 756 817 757 818 #else
Note: See TracChangeset
for help on using the changeset viewer.