- Timestamp:
- 2017-09-05T19:53:41+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90
r8496 r8498 157 157 !!------------------------------------------------------------------ 158 158 INTEGER :: ji, jj, jk, jl ! dummy loop indices 159 REAL(wp) :: ze_i, z1_ CpR, zdiscrim, zaaa, z1_2aaa! local scalars160 REAL(wp) :: ze_s, z L_Cp , ztmelts , zbbb, zccc! - -161 REAL(wp) :: zhmax, z1_zhmax, zsm_i , zcpMcpic, zt_i! - -162 REAL(wp) :: zlay_i, zlay_s ! - -159 REAL(wp) :: ze_i, z1_cp, z1_2cp ! local scalars 160 REAL(wp) :: ze_s, ztmelts, zbbb, zccc ! - - 161 REAL(wp) :: zhmax, z1_zhmax, zsm_i ! - - 162 REAL(wp) :: zlay_i, zlay_s ! - - 163 163 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i, z1_v_i 164 164 !!------------------------------------------------------------------ … … 181 181 ht_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:) !--- ice thickness 182 182 183 183 zhmax = hi_max(jpl) 184 184 z1_zhmax = 1._wp / hi_max(jpl) 185 185 WHERE( ht_i(:,:,jpl) > zhmax ) !--- bound ht_i by hi_max (i.e. 99 m) with associated update of ice area … … 205 205 !------------------- 206 206 zlay_i = REAL( nlay_i , wp ) ! number of layers 207 zaaa = cpic ! Conversion q(S,T) -> T (second order equation) 208 z1_2aaa = 1._wp / ( 2._wp * zaaa ) 209 zcpMcpic = rcp - cpic 207 z1_2cp = 1._wp / ( 2._wp * cpic ) 210 208 DO jl = 1, jpl 211 209 DO jk = 1, nlay_i … … 214 212 IF ( v_i(ji,jj,jl) > epsi20 ) THEN !--- icy area 215 213 ! 216 ze_i = e_i(ji,jj,jk,jl) / v_i(ji,jj,jl) * zlay_i ! Energy of melting e(S,T) [J.m-3] 217 ztmelts = - s_i(ji,jj,jk,jl) * tmut ! Ice layer melt temperature [C] 218 ! 219 zbbb = zcpMcpic * ztmelts + ze_i * r1_rhoic - lfus 220 zccc = lfus * ztmelts 221 zdiscrim = SQRT( MAX( zbbb*zbbb - 4._wp*zaaa*zccc , 0._wp) ) 222 zt_i = - ( zbbb + zdiscrim ) * z1_2aaa ! [C] 223 t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( zt_i , ztmelts ) ) + rt0 ! [K] with bounds: -100 < zt_i < ztmelts 214 ze_i = e_i(ji,jj,jk,jl) / v_i(ji,jj,jl) * zlay_i ! Energy of melting e(S,T) [J.m-3] 215 ztmelts = - s_i(ji,jj,jk,jl) * tmut ! Ice layer melt temperature [C] 216 ! Conversion q(S,T) -> T (second order equation) 217 zbbb = ( rcp - cpic ) * ztmelts + ze_i * r1_rhoic - lfus 218 zccc = SQRT( MAX( zbbb * zbbb - 4._wp * cpic * lfus * ztmelts , 0._wp) ) 219 t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * z1_2cp , ztmelts ) ) + rt0 ! [K] with bounds: -100 < t_i < ztmelts 224 220 ! 225 221 ELSE !--- no ice 226 t_i(ji,jj,jk,jl) = 222 t_i(ji,jj,jk,jl) = rt0 - 100._wp ! impose 173.15 K (i.e. -100 C) 227 223 !!clem: I think we should impose rt0 instead 228 224 ENDIF … … 235 231 ! Snow temperature [K] (with a minimum value (rt0 - 100.) imposed everywhere) 236 232 !-------------------- 237 z1_CpR = 1._wp / ( cpic * rhosn )238 zL_Cp = lfus / cpic239 233 zlay_s = REAL( nlay_s , wp ) 234 z1_cp = 1._wp / cpic 240 235 DO jk = 1, nlay_s 241 236 WHERE( v_s(:,:,:) > epsi20 ) !--- icy area 242 t_s(:,:,jk,:) = MAX( -100._wp , MIN( - z1_CpR * ( e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s ) + zL_Cp , 0._wp ) ) + rt0237 t_s(:,:,jk,:) = MAX( -100._wp , MIN( z1_cp * ( -r1_rhosn * (e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s) + lfus ) , 0._wp ) ) + rt0 243 238 ELSEWHERE !--- no ice 244 239 t_s(:,:,jk,:) = rt0 - 100._wp ! impose 173.15 K (i.e. -100 C) 245 240 END WHERE 246 241 END DO 247 !!gm perhaps better like this (?) : (if this compile, yes! one test instead of nlay_s tests)248 ! WHERE( v_s(:,:,:) > epsi20 ) !--- icy area249 ! DO jk = 1, nlay_s250 ! t_s(:,:,jk,:) = MAX( -100._wp , MIN( - z1_CpR * ( e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s ) + zL_Cp , 0._wp ) ) + rt0251 ! END DO252 ! ELSEWHERE !--- no ice253 ! DO jk = 1, nlay_s254 ! t_s(:,:,jk,:) = rt0 - 100._wp ! impose 173.15 K (i.e. -100 C)255 ! END DO256 ! END WHERE257 !!gm end better ?258 242 259 243 ! integrated values
Note: See TracChangeset
for help on using the changeset viewer.