Changeset 5260 for branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
- Timestamp:
- 2015-05-12T12:37:15+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r4333 r5260 19 19 USE phycst ! physical constants (ocean directory) 20 20 USE ice ! LIM-3 variables 21 USE par_ice ! LIM-3 parameters22 21 USE thd_ice ! LIM-3: thermodynamics 23 22 USE in_out_manager ! I/O manager … … 25 24 USE wrk_nemo ! work arrays 26 25 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 26 USE sbc_oce, ONLY : lk_cpl 27 27 28 28 IMPLICIT NONE … … 31 31 PUBLIC lim_thd_dif ! called by lim_thd 32 32 33 REAL(wp) :: epsi10 = 1.e-10_wp !34 33 !!---------------------------------------------------------------------- 35 34 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) … … 39 38 CONTAINS 40 39 41 SUBROUTINE lim_thd_dif( kideb , kiut , jl)40 SUBROUTINE lim_thd_dif( kideb , kiut ) 42 41 !!------------------------------------------------------------------ 43 42 !! *** ROUTINE lim_thd_dif *** … … 74 73 !! 75 74 !! ** Inputs / Ouputs : (global commons) 76 !! surface temperature : t_su_ b77 !! ice/snow temperatures : t_i_ b, t_s_b78 !! ice salinities : s_i_ b75 !! surface temperature : t_su_1d 76 !! ice/snow temperatures : t_i_1d, t_s_1d 77 !! ice salinities : s_i_1d 79 78 !! number of layers in the ice/snow: nlay_i, nlay_s 80 79 !! profile of the ice/snow layers : z_i, z_s 81 !! total ice/snow thickness : ht_i_ b, ht_s_b80 !! total ice/snow thickness : ht_i_1d, ht_s_1d 82 81 !! 83 82 !! ** External : … … 91 90 !! (04-2007) Energy conservation tested by M. Vancoppenolle 92 91 !!------------------------------------------------------------------ 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 92 INTEGER , INTENT(in) :: kideb, kiut ! Start/End point on which the the computation is applied 96 93 97 94 !! * Local variables … … 99 96 INTEGER :: ii, ij ! temporary dummy loop index 100 97 INTEGER :: numeq ! current reference number of equation 101 INTEGER :: layer! vertical dummy loop index98 INTEGER :: jk ! vertical dummy loop index 102 99 INTEGER :: nconv ! number of iterations in iterative procedure 103 100 INTEGER :: minnumeqmin, maxnumeqmax 104 INTEGER, DIMENSION(kiut) :: numeqmin ! reference number of top equation 105 INTEGER, DIMENSION(kiut) :: numeqmax ! reference number of bottom equation 106 INTEGER, DIMENSION(kiut) :: isnow ! switch for presence (1) or absence (0) of snow 101 102 INTEGER, POINTER, DIMENSION(:) :: numeqmin ! reference number of top equation 103 INTEGER, POINTER, DIMENSION(:) :: numeqmax ! reference number of bottom equation 104 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) :: zhsu 115 116 REAL(wp), POINTER, DIMENSION(:) :: isnow ! switch for presence (1) or absence (0) of snow 117 REAL(wp), POINTER, DIMENSION(:) :: ztsub ! old surface temperature (before the iterative procedure ) 118 REAL(wp), POINTER, DIMENSION(:) :: ztsubit ! surface temperature at previous iteration 119 REAL(wp), POINTER, DIMENSION(:) :: zh_i ! ice layer thickness 120 REAL(wp), POINTER, DIMENSION(:) :: zh_s ! snow layer thickness 121 REAL(wp), POINTER, DIMENSION(:) :: zfsw ! solar radiation absorbed at the surface 122 REAL(wp), POINTER, DIMENSION(:) :: zqns_ice_b ! solar radiation absorbed at the surface 123 REAL(wp), POINTER, DIMENSION(:) :: zf ! surface flux function 124 REAL(wp), POINTER, DIMENSION(:) :: dzf ! derivative of the surface flux function 125 REAL(wp), POINTER, DIMENSION(:) :: zerrit ! current error on temperature 126 REAL(wp), POINTER, DIMENSION(:) :: zdifcase ! case of the equation resolution (1->4) 127 REAL(wp), POINTER, DIMENSION(:) :: zftrice ! solar radiation transmitted through the ice 128 REAL(wp), POINTER, DIMENSION(:) :: zihic 129 130 REAL(wp), POINTER, DIMENSION(:,:) :: ztcond_i ! Ice thermal conductivity 131 REAL(wp), POINTER, DIMENSION(:,:) :: zradtr_i ! Radiation transmitted through the ice 132 REAL(wp), POINTER, DIMENSION(:,:) :: zradab_i ! Radiation absorbed in the ice 133 REAL(wp), POINTER, DIMENSION(:,:) :: zkappa_i ! Kappa factor in the ice 134 REAL(wp), POINTER, DIMENSION(:,:) :: ztib ! Old temperature in the ice 135 REAL(wp), POINTER, DIMENSION(:,:) :: zeta_i ! Eta factor in the ice 136 REAL(wp), POINTER, DIMENSION(:,:) :: ztitemp ! Temporary temperature in the ice to check the convergence 137 REAL(wp), POINTER, DIMENSION(:,:) :: zspeche_i ! Ice specific heat 138 REAL(wp), POINTER, DIMENSION(:,:) :: z_i ! Vertical cotes of the layers in the ice 139 REAL(wp), POINTER, DIMENSION(:,:) :: zradtr_s ! Radiation transmited through the snow 140 REAL(wp), POINTER, DIMENSION(:,:) :: zradab_s ! Radiation absorbed in the snow 141 REAL(wp), POINTER, DIMENSION(:,:) :: zkappa_s ! Kappa factor in the snow 142 REAL(wp), POINTER, DIMENSION(:,:) :: zeta_s ! Eta factor in the snow 143 REAL(wp), POINTER, DIMENSION(:,:) :: ztstemp ! Temporary temperature in the snow to check the convergence 144 REAL(wp), POINTER, DIMENSION(:,:) :: ztsb ! Temporary temperature in the snow 145 REAL(wp), POINTER, DIMENSION(:,:) :: z_s ! Vertical cotes of the layers in the snow 146 REAL(wp), POINTER, DIMENSION(:,:) :: zindterm ! 'Ind'ependent term 147 REAL(wp), POINTER, DIMENSION(:,:) :: zindtbis ! Temporary 'ind'ependent term 148 REAL(wp), POINTER, DIMENSION(:,:) :: zdiagbis ! Temporary 'dia'gonal term 149 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrid ! Tridiagonal system terms 150 151 ! diag errors on heat 152 REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini, zhfx_err 153 154 ! Mono-category 155 REAL(wp) :: zepsilon ! determines thres. above which computation of G(h) is done 156 REAL(wp) :: zratio_s ! dummy factor 157 REAL(wp) :: zratio_i ! dummy factor 158 REAL(wp) :: zh_thres ! thickness thres. for G(h) computation 159 REAL(wp) :: zhe ! dummy factor 160 REAL(wp) :: zkimean ! mean sea ice thermal conductivity 161 REAL(wp) :: zfac ! dummy factor 162 REAL(wp) :: zihe ! dummy factor 163 REAL(wp) :: zheshth ! dummy factor 164 165 REAL(wp), POINTER, DIMENSION(:) :: zghe ! G(he), th. conduct enhancement factor, mono-cat 166 147 167 !!------------------------------------------------------------------ 148 168 ! 169 CALL wrk_alloc( jpij, numeqmin, numeqmax ) 170 CALL wrk_alloc( jpij, isnow, ztsub, ztsubit, zh_i, zh_s, zfsw ) 171 CALL wrk_alloc( jpij, zf, dzf, zqns_ice_b, zerrit, zdifcase, zftrice, zihic, zghe ) 172 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 ) 173 CALL wrk_alloc( jpij,nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart=0 ) 174 CALL wrk_alloc( jpij,nlay_i+3, zindterm, zindtbis, zdiagbis ) 175 CALL wrk_alloc( jpij,nlay_i+3,3, ztrid ) 176 177 CALL wrk_alloc( jpij, zdq, zq_ini, zhfx_err ) 178 179 ! --- diag error on heat diffusion - PART 1 --- ! 180 zdq(:) = 0._wp ; zq_ini(:) = 0._wp 181 DO ji = kideb, kiut 182 zq_ini(ji) = ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i + & 183 & SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s ) 184 END DO 185 149 186 !------------------------------------------------------------------------------! 150 187 ! 1) Initialization ! 151 188 !------------------------------------------------------------------------------! 152 !153 189 DO ji = kideb , kiut 154 ! is there snow or not 155 isnow(ji)= NINT( 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_b(ji) ) ) ) 156 ! surface temperature of fusion 157 !!gm ??? ztfs(ji) = rtt !!!???? 158 ztfs(ji) = REAL( isnow(ji) ) * rtt + REAL( 1 - isnow(ji) ) * rtt 190 isnow(ji)= 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) ) ! is there snow or not 159 191 ! layer thickness 160 zh_i(ji) = ht_i_ b(ji) / REAL( nlay_i )161 zh_s(ji) = ht_s_ b(ji) / REAL( nlay_s )192 zh_i(ji) = ht_i_1d(ji) * r1_nlay_i 193 zh_s(ji) = ht_s_1d(ji) * r1_nlay_s 162 194 END DO 163 195 … … 169 201 z_i(:,0) = 0._wp ! vert. coord. of the up. lim. of the 1st ice layer 170 202 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 )203 DO jk = 1, nlay_s ! vert. coord of the up. lim. of the layer-th snow layer 204 DO ji = kideb , kiut 205 z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) * r1_nlay_s 206 END DO 207 END DO 208 209 DO jk = 1, nlay_i ! vert. coord of the up. lim. of the layer-th ice layer 210 DO ji = kideb , kiut 211 z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) * r1_nlay_i 180 212 END DO 181 213 END DO 182 214 ! 183 215 !------------------------------------------------------------------------------| 184 ! 2) Radiation s|216 ! 2) Radiation | 185 217 !------------------------------------------------------------------------------| 186 218 ! … … 194 226 ! zfsw = (1-i0).qsr_ice is absorbed at the surface 195 227 ! zftrice = io.qsr_ice is below the surface 196 ! fstbif = io.qsr_ice.exp(-k(h_i)) transmitted below the ice 197 228 ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice 229 ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover 230 zhsu = 0.1_wp ! threshold for the computation of i0 198 231 DO ji = kideb , kiut 199 232 ! switches 200 isnow(ji) = NINT( 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_b(ji) ) ))233 isnow(ji) = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) ) 201 234 ! hs > 0, isnow = 1 202 zhsu (ji) = hnzst ! threshold for the computation of i0 203 zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_b(ji) / zhsu(ji) ) ) 204 205 i0(ji) = REAL( 1 - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) ) 206 !fr1_i0_1d = i0 for a thin ice surface 207 !fr1_i0_2d = i0 for a thick ice surface 208 ! a function of the cloud cover 209 ! 210 !i0(ji) = (1.0-FLOAT(isnow(ji)))*3.0/(100*ht_s_b(ji)+10.0) 211 !formula used in Cice 235 zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_1d(ji) / zhsu ) ) 236 237 i0(ji) = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) ) 212 238 END DO 213 239 … … 217 243 !------------------------------------------------------- 218 244 DO ji = kideb , kiut 219 zfsw (ji) = qsr_ice_1d(ji) * ( 1 - i0(ji) ) ! Shortwave radiation absorbed at surface 220 zftrice(ji) = qsr_ice_1d(ji) * i0(ji) ! Solar radiation transmitted below the surface layer 221 dzf (ji) = dqns_ice_1d(ji) ! derivative of incoming nonsolar flux 245 zfsw (ji) = qsr_ice_1d(ji) * ( 1 - i0(ji) ) ! Shortwave radiation absorbed at surface 246 zftrice(ji) = qsr_ice_1d(ji) * i0(ji) ! Solar radiation transmitted below the surface layer 247 dzf (ji) = dqns_ice_1d(ji) ! derivative of incoming nonsolar flux 248 zqns_ice_b(ji) = qns_ice_1d(ji) ! store previous qns_ice_1d value 222 249 END DO 223 250 … … 230 257 END DO 231 258 232 DO layer= 1, nlay_s ! Radiation through snow259 DO jk = 1, nlay_s ! Radiation through snow 233 260 DO ji = kideb, kiut 234 261 ! ! 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) ) ) )262 zradtr_s(ji,jk) = zradtr_s(ji,0) * EXP( - zraext_s * ( MAX ( 0._wp , z_s(ji,jk) ) ) ) 236 263 ! ! radiation absorbed by the layer-th snow layer 237 zradab_s(ji, layer) = zradtr_s(ji,layer-1) - zradtr_s(ji,layer)264 zradab_s(ji,jk) = zradtr_s(ji,jk-1) - zradtr_s(ji,jk) 238 265 END DO 239 266 END DO 240 267 241 268 DO ji = kideb, kiut ! ice initialization 242 zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * REAL( isnow(ji) ) + zftrice(ji) * REAL( 1- isnow(ji) )243 END DO 244 245 DO layer= 1, nlay_i ! Radiation through ice269 zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * isnow(ji) + zftrice(ji) * ( 1._wp - isnow(ji) ) 270 END DO 271 272 DO jk = 1, nlay_i ! Radiation through ice 246 273 DO ji = kideb, kiut 247 274 ! ! 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) ) ) )275 zradtr_i(ji,jk) = zradtr_i(ji,0) * EXP( - rn_kappa_i * ( MAX ( 0._wp , z_i(ji,jk) ) ) ) 249 276 ! ! radiation absorbed by the layer-th ice layer 250 zradab_i(ji, layer) = zradtr_i(ji,layer-1) - zradtr_i(ji,layer)277 zradab_i(ji,jk) = zradtr_i(ji,jk-1) - zradtr_i(ji,jk) 251 278 END DO 252 279 END DO 253 280 254 281 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 282 ftr_ice_1d(ji) = zradtr_i(ji,nlay_i) 271 283 END DO 272 284 … … 277 289 ! 278 290 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 )! necessary282 zerrit (ji) = 1000._wp! initial value of error283 END DO 284 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)291 ztsub (ji) = t_su_1d(ji) ! temperature at the beg of iter pr. 292 ztsubit(ji) = t_su_1d(ji) ! temperature at the previous iter 293 t_su_1d(ji) = MIN( t_su_1d(ji), rt0 - ztsu_err ) ! necessary 294 zerrit (ji) = 1000._wp ! initial value of error 295 END DO 296 297 DO jk = 1, nlay_s ! Old snow temperature 298 DO ji = kideb , kiut 299 ztsb(ji,jk) = t_s_1d(ji,jk) 300 END DO 301 END DO 302 303 DO jk = 1, nlay_i ! Old ice temperature 304 DO ji = kideb , kiut 305 ztib(ji,jk) = t_i_1d(ji,jk) 294 306 END DO 295 307 END DO … … 298 310 zerritmax = 1000._wp ! maximal value of error on all points 299 311 300 DO WHILE ( zerritmax > maxer_i_thd .AND. nconv < nconv_i_thd)312 DO WHILE ( zerritmax > rn_terr_dif .AND. nconv < nn_conv_dif ) 301 313 ! 302 314 nconv = nconv + 1 … … 306 318 !------------------------------------------------------------------------------| 307 319 ! 308 IF( thcon_i_swi== 0 ) THEN ! Untersteiner (1964) formula309 DO ji = kideb , kiut 310 ztcond_i(ji,0) = rcdic + zbeta*s_i_b(ji,1) / MIN(-epsi10,t_i_b(ji,1)-rtt)311 ztcond_i(ji,0) = MAX(ztcond_i(ji,0),zkimin)312 END DO 313 DO layer= 1, nlay_i-1320 IF( nn_ice_thcon == 0 ) THEN ! Untersteiner (1964) formula 321 DO ji = kideb , kiut 322 ztcond_i(ji,0) = rcdic + zbeta * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) 323 ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin ) 324 END DO 325 DO jk = 1, nlay_i-1 314 326 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)327 ztcond_i(ji,jk) = rcdic + zbeta * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) / & 328 MIN(-2.0_wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0) 329 ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin ) 318 330 END DO 319 331 END DO 320 332 ENDIF 321 333 322 IF( thcon_i_swi== 1 ) THEN ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T323 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)334 IF( nn_ice_thcon == 1 ) THEN ! Pringle et al formula included: 2.11 + 0.09 S/T - 0.011.T 335 DO ji = kideb , kiut 336 ztcond_i(ji,0) = rcdic + 0.090_wp * s_i_1d(ji,1) / MIN( -epsi10, t_i_1d(ji,1) - rt0 ) & 337 & - 0.011_wp * ( t_i_1d(ji,1) - rt0 ) 326 338 ztcond_i(ji,0) = MAX( ztcond_i(ji,0), zkimin ) 327 339 END DO 328 DO layer= 1, nlay_i-1340 DO jk = 1, nlay_i-1 329 341 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 ) 342 ztcond_i(ji,jk) = rcdic + & 343 & 0.09_wp * ( s_i_1d(ji,jk) + s_i_1d(ji,jk+1) ) & 344 & / MIN( -2._wp * epsi10, t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0_wp * rt0 ) & 345 & - 0.0055_wp * ( t_i_1d(ji,jk) + t_i_1d(ji,jk+1) - 2.0 * rt0 ) 346 ztcond_i(ji,jk) = MAX( ztcond_i(ji,jk), zkimin ) 334 347 END DO 335 348 END DO 336 349 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)350 ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_1d(ji,nlay_i) / MIN( -epsi10, t_bo_1d(ji) - rt0 ) & 351 & - 0.011_wp * ( t_bo_1d(ji) - rt0 ) 339 352 ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin ) 340 353 END DO 341 354 ENDIF 342 ! 343 !------------------------------------------------------------------------------| 344 ! 5) kappa factors | 345 !------------------------------------------------------------------------------| 346 ! 355 356 ! 357 !------------------------------------------------------------------------------| 358 ! 5) G(he) - enhancement of thermal conductivity in mono-category case | 359 !------------------------------------------------------------------------------| 360 ! 361 ! Computation of effective thermal conductivity G(h) 362 ! Used in mono-category case only to simulate an ITD implicitly 363 ! Fichefet and Morales Maqueda, JGR 1997 364 365 zghe(:) = 1._wp 366 367 SELECT CASE ( nn_monocat ) 368 369 CASE (1,3) ! LIM3 370 371 zepsilon = 0.1_wp 372 zh_thres = EXP( 1._wp ) * zepsilon * 0.5_wp 373 374 DO ji = kideb, kiut 375 376 ! Mean sea ice thermal conductivity 377 zkimean = SUM( ztcond_i(ji,0:nlay_i) ) / REAL( nlay_i+1, wp ) 378 379 ! Effective thickness he (zhe) 380 zfac = 1._wp / ( rcdsn + zkimean ) 381 zratio_s = rcdsn * zfac 382 zratio_i = zkimean * zfac 383 zhe = zratio_s * ht_i_1d(ji) + zratio_i * ht_s_1d(ji) 384 385 ! G(he) 386 rswitch = MAX( 0._wp , SIGN( 1._wp , zhe - zh_thres ) ) ! =0 if zhe < zh_thres, if > 387 zghe(ji) = ( 1._wp - rswitch ) + rswitch * 0.5_wp * ( 1._wp + LOG( 2._wp * zhe / zepsilon ) ) 388 389 ! Impose G(he) < 2. 390 zghe(ji) = MIN( zghe(ji), 2._wp ) 391 392 END DO 393 394 END SELECT 395 396 ! 397 !------------------------------------------------------------------------------| 398 ! 6) kappa factors | 399 !------------------------------------------------------------------------------| 400 ! 401 !--- Snow 347 402 DO ji = kideb, kiut 348 349 !-- Snow kappa factors350 zkappa_s(ji, 0) = rcdsn / MAX(epsi10,zh_s(ji))351 zkappa_s(ji,nlay_s) = rcdsn / MAX(epsi10,zh_s(ji))352 END DO 353 354 DO layer = 1, nlay_s-1355 DO ji = kideb , kiut356 zkappa_s(ji,layer) = 2.0 * rcdsn / &357 MAX(epsi10,2.0*zh_s(ji))358 END DO 359 END DO360 361 DO layer = 1, nlay_i-1362 DO ji = kideb , kiut363 !-- Ice kappa factors364 zkappa_i(ji,layer) = 2.0*ztcond_i(ji,layer)/ &365 MAX(epsi10,2.0*zh_i(ji)) 366 END DO367 END DO368 369 DO ji = kideb , kiut370 zkappa_i(ji, 0) = ztcond_i(ji,0)/MAX(epsi10,zh_i(ji))371 zkappa_ i(ji,nlay_i) = ztcond_i(ji,nlay_i) / MAX(epsi10,zh_i(ji))372 !-- Interface373 zkappa_ s(ji,nlay_s) = 2.0*rcdsn*ztcond_i(ji,0)/MAX(epsi10, &374 (ztcond_i(ji,0)*zh_s(ji) + rcdsn*zh_i(ji)))375 zkappa_i(ji,0) = zkappa_s(ji,nlay_s)*REAL( isnow(ji) ) & 376 + zkappa_i(ji,0)*REAL( 1 - isnow(ji) )377 END DO378 ! 379 !------------------------------------------------------------------------------| 380 ! 6) Sea ice specific heat, eta factors |381 !------------------------------------------------------------------------------|382 !383 DO layer = 1, nlay_i384 DO ji = kideb , kiut385 z titemp(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), & 389 epsi10)390 END DO391 END DO392 393 DO layer = 1, nlay_s394 DO ji = kideb , kiut395 ztstemp(ji,layer) = t_s_b(ji,layer) 396 zeta_s(ji,layer) = rdt_ice / MAX(rhosn*cpic*zh_s(ji),epsi10)397 END DO398 END DO399 ! 400 ! ------------------------------------------------------------------------------|401 ! 7) surface flux computation |402 !------------------------------------------------------------------------------|403 !404 DO ji = kideb , kiut405 406 ! update of the non solar flux according to the update in T_su407 qnsr_ice_1d(ji) = qnsr_ice_1d(ji) + dqns_ice_1d(ji) * & 408 ( t_su_b(ji) - ztsuoldit(ji) )409 403 zfac = 1. / MAX( epsi10 , zh_s(ji) ) 404 zkappa_s(ji,0) = zghe(ji) * rcdsn * zfac 405 zkappa_s(ji,nlay_s) = zghe(ji) * rcdsn * zfac 406 END DO 407 408 DO jk = 1, nlay_s-1 409 DO ji = kideb , kiut 410 zkappa_s(ji,jk) = zghe(ji) * 2.0 * rcdsn / MAX( epsi10, 2.0 * zh_s(ji) ) 411 END DO 412 END DO 413 414 !--- Ice 415 DO jk = 1, nlay_i-1 416 DO ji = kideb , kiut 417 zkappa_i(ji,jk) = zghe(ji) * 2.0 * ztcond_i(ji,jk) / MAX( epsi10 , 2.0 * zh_i(ji) ) 418 END DO 419 END DO 420 421 !--- Snow-ice interface 422 DO ji = kideb , kiut 423 zfac = 1./ MAX( epsi10 , zh_i(ji) ) 424 zkappa_i(ji,0) = zghe(ji) * ztcond_i(ji,0) * zfac 425 zkappa_i(ji,nlay_i) = zghe(ji) * ztcond_i(ji,nlay_i) * zfac 426 zkappa_s(ji,nlay_s) = zghe(ji) * zghe(ji) * 2.0 * rcdsn * ztcond_i(ji,0) / & 427 & MAX( epsi10, ( zghe(ji) * ztcond_i(ji,0) * zh_s(ji) + zghe(ji) * rcdsn * zh_i(ji) ) ) 428 zkappa_i(ji,0) = zkappa_s(ji,nlay_s) * isnow(ji) + zkappa_i(ji,0) * ( 1._wp - isnow(ji) ) 429 END DO 430 431 ! 432 !------------------------------------------------------------------------------| 433 ! 7) Sea ice specific heat, eta factors | 434 !------------------------------------------------------------------------------| 435 ! 436 DO jk = 1, nlay_i 437 DO ji = kideb , kiut 438 ztitemp(ji,jk) = t_i_1d(ji,jk) 439 zspeche_i(ji,jk) = cpic + zgamma * s_i_1d(ji,jk) / MAX( ( t_i_1d(ji,jk) - rt0 ) * ( ztib(ji,jk) - rt0 ), epsi10 ) 440 zeta_i(ji,jk) = rdt_ice / MAX( rhoic * zspeche_i(ji,jk) * zh_i(ji), epsi10 ) 441 END DO 442 END DO 443 444 DO jk = 1, nlay_s 445 DO ji = kideb , kiut 446 ztstemp(ji,jk) = t_s_1d(ji,jk) 447 zeta_s(ji,jk) = rdt_ice / MAX( rhosn * cpic * zh_s(ji), epsi10 ) 448 END DO 449 END DO 450 451 ! 452 !------------------------------------------------------------------------------| 453 ! 8) surface flux computation | 454 !------------------------------------------------------------------------------| 455 ! 456 IF ( ln_it_qnsice ) THEN 457 DO ji = kideb , kiut 458 ! update of the non solar flux according to the update in T_su 459 qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsubit(ji) ) 460 END DO 461 ENDIF 462 463 ! Update incoming flux 464 DO ji = kideb , kiut 410 465 ! update incoming flux 411 zf(ji) = zfsw(ji) & ! net absorbed solar radiation 412 + qnsr_ice_1d(ji) ! non solar total flux 413 ! (LWup, LWdw, SH, LH) 414 415 END DO 416 417 ! 418 !------------------------------------------------------------------------------| 419 ! 8) tridiagonal system terms | 466 zf(ji) = zfsw(ji) & ! net absorbed solar radiation 467 & + qns_ice_1d(ji) ! non solar total flux (LWup, LWdw, SH, LH) 468 END DO 469 470 ! 471 !------------------------------------------------------------------------------| 472 ! 9) tridiagonal system terms | 420 473 !------------------------------------------------------------------------------| 421 474 ! … … 427 480 !!ice interior terms (top equation has the same form as the others) 428 481 429 DO numeq=1, jkmax+2482 DO numeq=1,nlay_i+3 430 483 DO ji = kideb , kiut 431 484 ztrid(ji,numeq,1) = 0. … … 440 493 DO numeq = nlay_s + 2, nlay_s + nlay_i 441 494 DO ji = kideb , kiut 442 layer = numeq - nlay_s - 1 443 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 zindterm(ji,numeq) = ztiold(ji,layer) + zeta_i(ji,layer)* & 448 zradab_i(ji,layer) 495 jk = numeq - nlay_s - 1 496 ztrid(ji,numeq,1) = - zeta_i(ji,jk) * zkappa_i(ji,jk-1) 497 ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) ) 498 ztrid(ji,numeq,3) = - zeta_i(ji,jk) * zkappa_i(ji,jk) 499 zindterm(ji,numeq) = ztib(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk) 449 500 END DO 450 501 ENDDO … … 454 505 !!ice bottom term 455 506 ztrid(ji,numeq,1) = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1) 456 ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,nlay_i)*( zkappa_i(ji,nlay_i)*zg1 & 457 + zkappa_i(ji,nlay_i-1) ) 507 ztrid(ji,numeq,2) = 1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) ) 458 508 ztrid(ji,numeq,3) = 0.0 459 zindterm(ji,numeq) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i)* & 460 ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i)*zg1 & 461 * t_bo_b(ji) ) 509 zindterm(ji,numeq) = ztib(ji,nlay_i) + zeta_i(ji,nlay_i) * & 510 & ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 * t_bo_1d(ji) ) 462 511 ENDDO 463 512 464 513 465 514 DO ji = kideb , kiut 466 IF ( ht_s_ b(ji).gt.0.0 ) THEN515 IF ( ht_s_1d(ji) > 0.0 ) THEN 467 516 ! 468 517 !------------------------------------------------------------------------------| … … 472 521 !!snow interior terms (bottom equation has the same form as the others) 473 522 DO numeq = 3, nlay_s + 1 474 layer = numeq - 1 475 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 zindterm(ji,numeq) = ztsold(ji,layer) + zeta_s(ji,layer)* & 480 zradab_s(ji,layer) 523 jk = numeq - 1 524 ztrid(ji,numeq,1) = - zeta_s(ji,jk) * zkappa_s(ji,jk-1) 525 ztrid(ji,numeq,2) = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) ) 526 ztrid(ji,numeq,3) = - zeta_s(ji,jk)*zkappa_s(ji,jk) 527 zindterm(ji,numeq) = ztsb(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 481 528 END DO 482 529 … … 484 531 IF ( nlay_i.eq.1 ) THEN 485 532 ztrid(ji,nlay_s+2,3) = 0.0 486 zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1)* & 487 t_bo_b(ji) 533 zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji) 488 534 ENDIF 489 535 490 IF ( t_su_ b(ji) .LT. rtt) THEN536 IF ( t_su_1d(ji) < rt0 ) THEN 491 537 492 538 !------------------------------------------------------------------------------| … … 498 544 499 545 !!surface equation 500 ztrid(ji,1,1) = 0.0501 ztrid(ji,1,2) = dzf(ji) - zg1s*zkappa_s(ji,0)502 ztrid(ji,1,3) = zg1s*zkappa_s(ji,0)503 zindterm(ji,1) = dzf(ji) *t_su_b(ji)- zf(ji)546 ztrid(ji,1,1) = 0.0 547 ztrid(ji,1,2) = dzf(ji) - zg1s * zkappa_s(ji,0) 548 ztrid(ji,1,3) = zg1s * zkappa_s(ji,0) 549 zindterm(ji,1) = dzf(ji) * t_su_1d(ji) - zf(ji) 504 550 505 551 !!first layer of snow equation 506 ztrid(ji,2,1) = - zkappa_s(ji,0) *zg1s*zeta_s(ji,1)507 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) *(zkappa_s(ji,1) + zkappa_s(ji,0)*zg1s)552 ztrid(ji,2,1) = - zkappa_s(ji,0) * zg1s * zeta_s(ji,1) 553 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 508 554 ztrid(ji,2,3) = - zeta_s(ji,1)* zkappa_s(ji,1) 509 zindterm(ji,2) = zts old(ji,1) + zeta_s(ji,1)*zradab_s(ji,1)555 zindterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) * zradab_s(ji,1) 510 556 511 557 ELSE … … 521 567 !!first layer of snow equation 522 568 ztrid(ji,2,1) = 0.0 523 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + & 524 zkappa_s(ji,0) * zg1s ) 569 ztrid(ji,2,2) = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 525 570 ztrid(ji,2,3) = - zeta_s(ji,1)*zkappa_s(ji,1) 526 zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * & 527 ( zradab_s(ji,1) + & 528 zkappa_s(ji,0) * zg1s * t_su_b(ji) ) 571 zindterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) * & 572 & ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 529 573 ENDIF 530 574 ELSE … … 534 578 !------------------------------------------------------------------------------| 535 579 ! 536 IF ( t_su_b(ji) .LT. rtt) THEN580 IF ( t_su_1d(ji) < rt0 ) THEN 537 581 ! 538 582 !------------------------------------------------------------------------------| … … 548 592 ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0)*zg1 549 593 ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0)*zg1 550 zindterm(ji,numeqmin(ji)) = dzf(ji)*t_su_ b(ji) - zf(ji)594 zindterm(ji,numeqmin(ji)) = dzf(ji)*t_su_1d(ji) - zf(ji) 551 595 552 596 !!first layer of ice equation 553 597 ztrid(ji,numeqmin(ji)+1,1) = - zkappa_i(ji,0) * zg1 * zeta_i(ji,1) 554 ztrid(ji,numeqmin(ji)+1,2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) & 555 + zkappa_i(ji,0) * zg1 ) 556 ztrid(ji,numeqmin(ji)+1,3) = - zeta_i(ji,1)*zkappa_i(ji,1) 557 zindterm(ji,numeqmin(ji)+1)= ztiold(ji,1) + zeta_i(ji,1)*zradab_i(ji,1) 598 ztrid(ji,numeqmin(ji)+1,2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 599 ztrid(ji,numeqmin(ji)+1,3) = - zeta_i(ji,1) * zkappa_i(ji,1) 600 zindterm(ji,numeqmin(ji)+1)= ztib(ji,1) + zeta_i(ji,1) * zradab_i(ji,1) 558 601 559 602 !!case of only one layer in the ice (surface & ice equations are altered) 560 603 561 IF ( nlay_i.eq.1) THEN604 IF ( nlay_i == 1 ) THEN 562 605 ztrid(ji,numeqmin(ji),1) = 0.0 563 ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0)*2.0 564 ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0)*2.0 565 ztrid(ji,numeqmin(ji)+1,1) = -zkappa_i(ji,0)*2.0*zeta_i(ji,1) 566 ztrid(ji,numeqmin(ji)+1,2) = 1.0 + zeta_i(ji,1)*(zkappa_i(ji,0)*2.0 + & 567 zkappa_i(ji,1)) 606 ztrid(ji,numeqmin(ji),2) = dzf(ji) - zkappa_i(ji,0) * 2.0 607 ztrid(ji,numeqmin(ji),3) = zkappa_i(ji,0) * 2.0 608 ztrid(ji,numeqmin(ji)+1,1) = -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1) 609 ztrid(ji,numeqmin(ji)+1,2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) 568 610 ztrid(ji,numeqmin(ji)+1,3) = 0.0 569 611 570 zindterm(ji,numeqmin(ji)+1) = zti old(ji,1) + zeta_i(ji,1)*&571 ( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_b(ji) )612 zindterm(ji,numeqmin(ji)+1) = ztib(ji,1) + zeta_i(ji,1) * & 613 & ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) 572 614 ENDIF 573 615 … … 585 627 !!first layer of ice equation 586 628 ztrid(ji,numeqmin(ji),1) = 0.0 587 ztrid(ji,numeqmin(ji),2) = 1.0 + zeta_i(ji,1)*(zkappa_i(ji,1) + zkappa_i(ji,0)* & 588 zg1) 629 ztrid(ji,numeqmin(ji),2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 589 630 ztrid(ji,numeqmin(ji),3) = - zeta_i(ji,1) * zkappa_i(ji,1) 590 zindterm(ji,numeqmin(ji)) = zti old(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) +&591 zkappa_i(ji,0) * zg1 * t_su_b(ji) )631 zindterm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1) * & 632 & ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 592 633 593 634 !!case of only one layer in the ice (surface & ice equations are altered) 594 IF ( nlay_i.eq.1) THEN635 IF ( nlay_i == 1 ) THEN 595 636 ztrid(ji,numeqmin(ji),1) = 0.0 596 ztrid(ji,numeqmin(ji),2) = 1.0 + zeta_i(ji,1)*(zkappa_i(ji,0)*2.0 + & 597 zkappa_i(ji,1)) 637 ztrid(ji,numeqmin(ji),2) = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) 598 638 ztrid(ji,numeqmin(ji),3) = 0.0 599 zindterm(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.0 639 zindterm(ji,numeqmin(ji)) = ztib(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) & 640 & + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0 602 641 ENDIF 603 642 … … 609 648 ! 610 649 !------------------------------------------------------------------------------| 611 ! 9) tridiagonal system solving|650 ! 10) tridiagonal system solving | 612 651 !------------------------------------------------------------------------------| 613 652 ! … … 618 657 619 658 maxnumeqmax = 0 620 minnumeqmin = jkmax+4659 minnumeqmin = nlay_i+5 621 660 622 661 DO ji = kideb , kiut … … 627 666 END DO 628 667 629 DO layer = minnumeqmin+1, maxnumeqmax 630 DO ji = kideb , kiut 631 numeq = min(max(numeqmin(ji)+1,layer),numeqmax(ji)) 632 zdiagbis(ji,numeq) = ztrid(ji,numeq,2) - ztrid(ji,numeq,1)* & 633 ztrid(ji,numeq-1,3)/zdiagbis(ji,numeq-1) 634 zindtbis(ji,numeq) = zindterm(ji,numeq) - ztrid(ji,numeq,1)* & 635 zindtbis(ji,numeq-1)/zdiagbis(ji,numeq-1) 668 DO jk = minnumeqmin+1, maxnumeqmax 669 DO ji = kideb , kiut 670 numeq = min(max(numeqmin(ji)+1,jk),numeqmax(ji)) 671 zdiagbis(ji,numeq) = ztrid(ji,numeq,2) - ztrid(ji,numeq,1) * ztrid(ji,numeq-1,3) / zdiagbis(ji,numeq-1) 672 zindtbis(ji,numeq) = zindterm(ji,numeq) - ztrid(ji,numeq,1) * zindtbis(ji,numeq-1) / zdiagbis(ji,numeq-1) 636 673 END DO 637 674 END DO … … 639 676 DO ji = kideb , kiut 640 677 ! ice temperatures 641 t_i_b(ji,nlay_i) = zindtbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji)) 642 END DO 643 644 DO numeq = nlay_i + nlay_s + 1, nlay_s + 2, -1 645 DO ji = kideb , kiut 646 layer = numeq - nlay_s - 1 647 t_i_b(ji,layer) = (zindtbis(ji,numeq) - ztrid(ji,numeq,3)* & 648 t_i_b(ji,layer+1))/zdiagbis(ji,numeq) 678 t_i_1d(ji,nlay_i) = zindtbis(ji,numeqmax(ji)) / zdiagbis(ji,numeqmax(ji)) 679 END DO 680 681 DO numeq = nlay_i + nlay_s, nlay_s + 2, -1 682 DO ji = kideb , kiut 683 jk = numeq - nlay_s - 1 684 t_i_1d(ji,jk) = ( zindtbis(ji,numeq) - ztrid(ji,numeq,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,numeq) 649 685 END DO 650 686 END DO … … 652 688 DO ji = kideb , kiut 653 689 ! 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))) 690 IF (ht_s_1d(ji) > 0._wp) & 691 t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) & 692 & / zdiagbis(ji,nlay_s+1) * MAX( 0.0, SIGN( 1.0, ht_s_1d(ji) ) ) 658 693 659 694 ! 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))695 isnow(ji) = 1._wp - MAX( 0._wp , SIGN( 1._wp , -ht_s_1d(ji) ) ) 696 ztsubit(ji) = t_su_1d(ji) 697 IF( t_su_1d(ji) < rt0 ) & 698 t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3) * & 699 & ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji)) 665 700 END DO 666 701 ! 667 702 !-------------------------------------------------------------------------- 668 ! 1 0) Has the scheme converged ?, end of the iterative procedure |703 ! 11) Has the scheme converged ?, end of the iterative procedure | 669 704 !-------------------------------------------------------------------------- 670 705 ! 671 706 ! check that nowhere it has started to melt 672 ! zerrit(ji) is a measure of error, it has to be under maxer_i_thd 673 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))) 707 ! zerrit(ji) is a measure of error, it has to be under terr_dif 708 DO ji = kideb , kiut 709 t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , 190._wp ) 710 zerrit(ji) = ABS( t_su_1d(ji) - ztsubit(ji) ) 711 END DO 712 713 DO jk = 1, nlay_s 714 DO ji = kideb , kiut 715 t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), 190._wp ) 716 zerrit(ji) = MAX( zerrit(ji), ABS( t_s_1d(ji,jk) - ztstemp(ji,jk) ) ) 717 END DO 718 END DO 719 720 DO jk = 1, nlay_i 721 DO ji = kideb , kiut 722 ztmelt_i = -tmut * s_i_1d(ji,jk) + rt0 723 t_i_1d(ji,jk) = MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), 190._wp ) 724 zerrit(ji) = MAX( zerrit(ji), ABS( t_i_1d(ji,jk) - ztitemp(ji,jk) ) ) 692 725 END DO 693 726 END DO … … 703 736 END DO ! End of the do while iterative procedure 704 737 705 IF( ln_ nicep.AND. lwp ) THEN738 IF( ln_icectl .AND. lwp ) THEN 706 739 WRITE(numout,*) ' zerritmax : ', zerritmax 707 740 WRITE(numout,*) ' nconv : ', nconv … … 710 743 ! 711 744 !-------------------------------------------------------------------------! 712 ! 1 1) Fluxes at the interfaces !745 ! 12) Fluxes at the interfaces ! 713 746 !-------------------------------------------------------------------------! 714 747 DO ji = kideb, kiut 715 #if ! defined key_coupled716 748 ! 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 749 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 750 ! ! 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))751 isnow(ji) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_1d(ji) ) ) 752 fc_su(ji) = - isnow(ji) * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji)) & 753 & - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1 * (t_i_1d(ji,1) - t_su_1d(ji)) 723 754 ! ! 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 755 fc_bo_i(ji) = - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) ) 756 END DO 757 758 ! --- computes sea ice energy of melting compulsory for limthd_dh --- ! 759 CALL lim_thd_enmelt( kideb, kiut ) 760 761 ! --- diagnose the change in non-solar flux due to surface temperature change --- ! 762 IF ( ln_it_qnsice ) hfx_err_dif_1d(:) = hfx_err_dif_1d(:) - ( qns_ice_1d(:) - zqns_ice_b(:) ) * a_i_1d(:) 763 764 ! --- diag conservation imbalance on heat diffusion - PART 2 --- ! 765 DO ji = kideb, kiut 766 zdq(ji) = - zq_ini(ji) + ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i + & 767 & SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s ) 768 IF( t_su_1d(ji) < rt0 ) THEN ! case T_su < 0degC 769 zhfx_err(ji) = qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice 770 ELSE ! case T_su = 0degC 771 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 772 ENDIF 773 hfx_err_1d(ji) = hfx_err_1d(ji) + zhfx_err(ji) * a_i_1d(ji) 774 775 ! total heat that is sent to the ocean (i.e. not used in the heat diffusion equation) 776 hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err(ji) * a_i_1d(ji) 777 END DO 778 779 !----------------------------------------- 780 ! Heat flux used to warm/cool ice in W.m-2 781 !----------------------------------------- 782 DO ji = kideb, kiut 783 IF( t_su_1d(ji) < rt0 ) THEN ! case T_su < 0degC 784 hfx_dif_1d(ji) = hfx_dif_1d(ji) + & 785 & ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji) 786 ELSE ! case T_su = 0degC 787 hfx_dif_1d(ji) = hfx_dif_1d(ji) + & 788 & ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_1d(ji) 789 ENDIF 790 ! correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation 791 hfx_dif_1d(ji) = hfx_dif_1d(ji) - zhfx_err(ji) * a_i_1d(ji) 792 END DO 793 ! 794 CALL wrk_dealloc( jpij, numeqmin, numeqmax ) 795 CALL wrk_dealloc( jpij, isnow, ztsub, ztsubit, zh_i, zh_s, zfsw ) 796 CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zghe ) 797 CALL wrk_dealloc( jpij,nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 ) 798 CALL wrk_dealloc( jpij,nlay_s+1, zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart = 0 ) 799 CALL wrk_dealloc( jpij,nlay_i+3, zindterm, zindtbis, zdiagbis ) 800 CALL wrk_dealloc( jpij,nlay_i+3,3, ztrid ) 801 CALL wrk_dealloc( jpij, zdq, zq_ini, zhfx_err ) 802 803 END SUBROUTINE lim_thd_dif 804 805 SUBROUTINE lim_thd_enmelt( kideb, kiut ) 806 !!----------------------------------------------------------------------- 807 !! *** ROUTINE lim_thd_enmelt *** 808 !! 809 !! ** Purpose : Computes sea ice energy of melting q_i (J.m-3) from temperature 810 !! 811 !! ** Method : Formula (Bitz and Lipscomb, 1999) 812 !!------------------------------------------------------------------- 813 INTEGER, INTENT(in) :: kideb, kiut ! bounds for the spatial loop 814 ! 815 INTEGER :: ji, jk ! dummy loop indices 816 REAL(wp) :: ztmelts ! local scalar 817 !!------------------------------------------------------------------- 818 ! 819 DO jk = 1, nlay_i ! Sea ice energy of melting 731 820 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 821 ztmelts = - tmut * s_i_1d(ji,jk) + rt0 822 t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts ) ! Force t_i_1d to be lower than melting point 823 ! (sometimes dif scheme produces abnormally high temperatures) 824 q_i_1d(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_1d(ji,jk) ) & 825 & + lfus * ( 1.0 - ( ztmelts-rt0 ) / ( t_i_1d(ji,jk) - rt0 ) ) & 826 & - rcp * ( ztmelts-rt0 ) ) 827 END DO 828 END DO 829 DO jk = 1, nlay_s ! Snow energy of melting 830 DO ji = kideb, kiut 831 q_s_1d(ji,jk) = rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) 832 END DO 833 END DO 754 834 ! 755 END SUBROUTINE lim_thd_ dif835 END SUBROUTINE lim_thd_enmelt 756 836 757 837 #else
Note: See TracChangeset
for help on using the changeset viewer.
