- 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/limupdate2.F90
r4332 r4634 45 45 PUBLIC lim_update2 ! routine called by ice_step 46 46 47 REAL(wp) :: epsi10 = 1.e-10_wp ! - - 48 REAL(wp) :: rzero = 0._wp ! - - 49 REAL(wp) :: rone = 1._wp ! - - 50 47 REAL(wp) :: epsi10 = 1.e-10_wp ! - - 48 REAL(wp) :: epsi20 = 1.e-20_wp 49 51 50 !! * Substitutions 52 51 # include "vectopt_loop_substitute.h90" … … 77 76 INTEGER :: jbnd1, jbnd2 78 77 INTEGER :: i_ice_switch 79 INTEGER :: ind_im, layer ! indices for internal melt80 REAL(wp) :: z weight, zesum, zhimax, z_da_i81 REAL(wp) :: zinda, zindb, zindsn, zindic 82 REAL(wp) :: z indg, zh, zdvres, zviold283 REAL(wp) :: z bigvalue, zvsold2, z_da_ex84 REAL(wp) :: z _prescr_hi, zat_i_old, ztmelts, ze_s85 86 INTEGER , POINTER, DIMENSION(:,:,:) :: internal_melt87 REAL(wp), POINTER, DIMENSION(:) :: zthick0, zqm0 ! thickness of the layers and heat contents for 88 REAL(wp) :: zchk_v_i, zchk_smv, zchk_ fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset)78 REAL(wp) :: zindb, zindsn, zindic 79 REAL(wp) :: zh, zdvres, zsal 80 81 REAL(wp) :: zEs ! specific enthalpy of snow (J/kg) 82 REAL(wp) :: zEi ! specific enthalpy of ice (J/kg) 83 REAL(wp) :: zEw ! specific enthalpy of exchanged water (J/kg) 84 REAL(wp) :: zdE ! specific enthalpy difference (J/kg) 85 REAL(wp) :: zfmdt ! exchange mass flux x time step (J/m2), >0 towards the ocean 86 87 REAL(wp) :: zchk_v_i, zchk_smv, zchk_e_i, zchk_fs, zchk_fw, zchk_ft, zchk_v_i_b, zchk_smv_b, zchk_e_i_b, zchk_fs_b, zchk_fw_b, zchk_ft_b ! Check conservation (C Rousset) 89 88 REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 90 ! mass and salt flux (clem)91 REAL(wp), POINTER, DIMENSION(:,:,:) :: zviold, zvsold, zsmvold ! old ice volume...92 89 !!------------------------------------------------------------------- 93 90 IF( nn_timing == 1 ) CALL timing_start('limupdate2') 94 95 CALL wrk_alloc( jpi,jpj,jpl, internal_melt ) ! integer96 CALL wrk_alloc( jkmax, zthick0, zqm0 )97 98 CALL wrk_alloc( jpi,jpj,jpl,zviold, zvsold, zsmvold ) ! clem99 91 100 92 !---------------------------------------------------------------------------------------- 101 93 ! 1. Computation of trend terms 102 94 !---------------------------------------------------------------------------------------- 103 !- Trend terms104 d_a_i_thd(:,:,:) = a_i(:,:,:) - old_a_i(:,:,:)105 d_v_s_thd(:,:,:) = v_s(:,:,:) - old_v_s(:,:,:)106 d_v_i_thd(:,:,:) = v_i(:,:,:) - old_v_i(:,:,:)107 d_e_s_thd(:,:,:,:) = e_s(:,:,:,:) - old_e_s(:,:,:,:)108 d_e_i_thd(:,:,:,:) = e_i(:,:,:,:) - old_e_i(:,:,:,:)109 !?? d_oa_i_thd(:,:,:) = oa_i (:,:,:) - old_oa_i (:,:,:)110 d_smv_i_thd(:,:,:) = 0._wp111 IF( num_sal == 2 ) d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:)112 ! diag only (clem)113 dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday114 115 ! mass and salt flux init (clem)116 zviold(:,:,:) = v_i(:,:,:)117 zvsold(:,:,:) = v_s(:,:,:)118 zsmvold(:,:,:) = smv_i(:,:,:)119 95 120 96 ! ------------------------------- 121 97 !- check conservation (C Rousset) 122 98 IF (ln_limdiahsb) THEN 123 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:) , dim=3 ) * area(:,:) * tms(:,:) )99 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) 124 100 zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 125 zchk_fw_b = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 126 zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 101 zchk_e_i_b = glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) 102 zchk_fw_b = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) 103 zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) 104 zchk_ft_b = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) 127 105 ENDIF 128 106 !- check conservation (C Rousset) 129 107 ! ------------------------------- 130 108 109 ! zap small values 110 !----------------- 111 CALL lim_itd_me_zapsmall 112 131 113 CALL lim_var_glo2eqv 132 114 … … 134 116 ! 2. Review of all pathological cases 135 117 !-------------------------------------- 136 137 ! clem: useless now138 !-------------------------------------------139 ! 2.1) Advection of ice in an ice-free cell140 !-------------------------------------------141 ! should be removed since it is treated after dynamics now142 ! zhimax = 5._wp143 ! ! first category144 ! DO jj = 1, jpj145 ! DO ji = 1, jpi146 ! !--- the thickness of such an ice is often out of bounds147 ! !--- thus we recompute a new area while conserving ice volume148 ! zat_i_old = SUM( old_a_i(ji,jj,:) )149 ! zindb = MAX( 0._wp, SIGN( 1._wp, ABS( d_a_i_thd(ji,jj,1) ) - epsi10 ) )150 ! IF ( ( ABS( d_v_i_thd(ji,jj,1) ) / MAX( ABS( d_a_i_thd(ji,jj,1) ),epsi10 ) * zindb .GT. zhimax ) &151 ! & .AND. ( ( v_i(ji,jj,1) / MAX( a_i(ji,jj,1), epsi10 ) * zindb ) .GT. zhimax ) &152 ! & .AND. ( zat_i_old .LT. 1.e-6 ) ) THEN ! new line153 ! ht_i(ji,jj,1) = hi_max(1) * 0.5_wp154 ! a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1)155 ! ENDIF156 ! END DO157 ! END DO158 159 ! zhimax = 20._wp160 ! ! other categories161 ! DO jl = 2, jpl162 ! jm = ice_types(jl)163 ! DO jj = 1, jpj164 ! DO ji = 1, jpi165 ! zindb = MAX( rzero, SIGN( rone, ABS( d_a_i_thd(ji,jj,jl)) - epsi10 ) )166 ! ! this correction is very tricky... sometimes, advection gets wrong i don't know why167 ! ! it makes problems when the advected volume and concentration do not seem to be168 ! ! related with each other169 ! ! the new thickness is sometimes very big!170 ! ! and sometimes d_a_i_trp and d_v_i_trp have different sign171 ! ! which of course is plausible172 ! ! but fuck! it fucks everything up :)173 ! IF ( ( ABS( d_v_i_thd(ji,jj,jl) ) / MAX( ABS( d_a_i_thd(ji,jj,jl) ), epsi10 ) * zindb .GT. zhimax ) &174 ! & .AND. ( v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zindb ) .GT. zhimax ) THEN175 ! ht_i(ji,jj,jl) = ( hi_max_typ(jl-ice_cat_bounds(jm,1),jm) + hi_max_typ(jl-ice_cat_bounds(jm,1)+1,jm) ) * 0.5_wp176 ! a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl)177 ! ENDIF178 ! END DO ! ji179 ! END DO !jj180 ! END DO !jl181 182 118 at_i(:,:) = 0._wp 183 119 DO jl = 1, jpl … … 194 130 END DO 195 131 196 !--------------------------------- 197 ! 2.3) Melt of an internal layer 198 !--------------------------------- 199 internal_melt(:,:,:) = 0 200 201 DO jl = 1, jpl 202 DO jk = 1, nlay_i 203 DO jj = 1, jpj 204 DO ji = 1, jpi 205 ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt 206 IF ( ( ( e_i(ji,jj,jk,jl) .LE. 0.0 ) .OR. ( t_i(ji,jj,jk,jl) .GE. ztmelts ) ) & 207 & .AND. ( v_i(ji,jj,jl) .GT. 0.0 ) .AND. ( a_i(ji,jj,jl) .GT. 0.0 ) ) THEN 208 internal_melt(ji,jj,jl) = 1 209 ENDIF 210 END DO ! ji 211 END DO ! jj 212 END DO !jk 213 END DO !jl 214 215 DO jl = 1, jpl 216 DO jj = 1, jpj 217 DO ji = 1, jpi 218 IF( internal_melt(ji,jj,jl) == 1 ) THEN 219 ! initial ice thickness 220 !----------------------- 221 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 222 223 ! reduce ice thickness 224 !----------------------- 225 ind_im = 0 226 zesum = 0.0 227 DO jk = 1, nlay_i 228 ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt 229 IF ( ( e_i(ji,jj,jk,jl) .LE. 0.0 ) .OR. ( t_i(ji,jj,jk,jl) .GE. ztmelts ) ) ind_im = ind_im + 1 230 zesum = zesum + e_i(ji,jj,jk,jl) 231 END DO 232 ht_i(ji,jj,jl) = ht_i(ji,jj,jl) - REAL(ind_im)*ht_i(ji,jj,jl) / REAL(nlay_i) 233 v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) 234 235 !CLEM 236 zdvres = REAL(ind_im)*ht_i(ji,jj,jl) / REAL(nlay_i) * a_i(ji,jj,jl) 237 !rdm_ice(ji,jj) = rdm_ice(ji,jj) - zdvres * rhoic 238 !sfx_res(ji,jj) = sfx_res(ji,jj) + sm_i(ji,jj,jl) * ( rhoic * zdvres / rdt_ice ) 239 240 ! redistribute heat 241 !----------------------- 242 ! old thicknesses and enthalpies 243 ind_im = 0 244 DO jk = 1, nlay_i 245 ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt 246 IF ( ( e_i(ji,jj,jk,jl) .GT. 0.0 ) .AND. & 247 ( t_i(ji,jj,jk,jl) .LT. ztmelts ) ) THEN 248 ind_im = ind_im + 1 249 zthick0(ind_im) = ht_i(ji,jj,jl) * REAL(ind_im / nlay_i) 250 zqm0 (ind_im) = MAX( e_i(ji,jj,jk,jl) , 0.0 ) 251 ENDIF 252 END DO 253 254 ! Redistributing energy on the new grid 255 IF ( ind_im .GT. 0 ) THEN 256 257 DO jk = 1, nlay_i 258 e_i(ji,jj,jk,jl) = 0.0 259 DO layer = 1, ind_im 260 zweight = MAX ( & 261 MIN( ht_i(ji,jj,jl) * REAL(layer/ind_im) , ht_i(ji,jj,jl) * REAL(jk / nlay_i) ) - & 262 MAX( ht_i(ji,jj,jl) * REAL((layer-1)/ind_im) , ht_i(ji,jj,jl) * REAL((jk-1) / nlay_i) ) , 0.0 ) & 263 / ( ht_i(ji,jj,jl) / REAL(ind_im) ) 264 265 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) + zweight*zqm0(layer) 266 END DO !layer 267 END DO ! jk 268 269 zesum = 0.0 270 DO jk = 1, nlay_i 271 zesum = zesum + e_i(ji,jj,jk,jl) 272 END DO 273 274 ELSE ! ind_im .EQ. 0, total melt 275 e_i(ji,jj,jk,jl) = 0.0 276 ENDIF 277 278 ENDIF ! internal_melt 279 280 END DO ! ji 281 END DO !jj 282 END DO !jl 283 284 internal_melt(:,:,:) = 0 285 286 287 ! Melt of snow 288 !-------------- 289 DO jl = 1, jpl 290 DO jj = 1, jpj 291 DO ji = 1, jpi 292 ! snow energy of melting 293 zinda = MAX( 0._wp, SIGN( 1._wp, v_s(ji,jj,jl) - epsi10 ) ) 294 ze_s = zinda * e_s(ji,jj,1,jl) * unit_fac / area(ji,jj) / MAX( v_s(ji,jj,jl), epsi10 ) ! snow energy of melting 295 296 ! If snow energy of melting smaller then Lf 297 ! Then all snow melts and meltwater, heat go to the ocean 298 IF ( ze_s .LE. rhosn * lfus ) internal_melt(ji,jj,jl) = 1 299 300 END DO 301 END DO 302 END DO 303 304 DO jl = 1, jpl 305 DO jj = 1, jpj 306 DO ji = 1, jpi 307 IF ( internal_melt(ji,jj,jl) == 1 ) THEN 308 zdvres = v_s(ji,jj,jl) 309 ! release heat 310 fheat_res(ji,jj) = fheat_res(ji,jj) + ze_s * zdvres / rdt_ice 311 ! release mass 312 !rdm_snw(ji,jj) = rdm_snw(ji,jj) - zdvres * rhosn 313 ! 314 v_s(ji,jj,jl) = 0.0 315 e_s(ji,jj,1,jl) = 0.0 316 ENDIF 317 END DO 318 END DO 319 END DO 320 321 zbigvalue = 1.0e+20 322 DO jl = 1, jpl 323 DO jj = 1, jpj 324 DO ji = 1, jpi 325 326 !switches 327 zindb = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi10 ) ) 328 !switch = 1 if a_i > 1e-06 and 0 if not 329 zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) ) !=1 if hs > 1e-10 and 0 if not 330 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) !=1 if hi > 1e-10 and 0 if not 331 ! bug fix 25 avril 2007 332 zindb = zindb*zindic 333 334 !--- 2.3 Correction to ice age 335 !------------------------------ 336 ! IF ((o_i(ji,jj,jl)-1.0)*rday.gt.(rdt_ice*float(numit))) THEN 337 ! o_i(ji,jj,jl) = rdt_ice*FLOAT(numit)/rday 338 ! ENDIF 339 IF ((oa_i(ji,jj,jl)-1.0)*rday.gt.(rdt_ice*numit*a_i(ji,jj,jl))) THEN 340 oa_i(ji,jj,jl) = rdt_ice*numit/rday*a_i(ji,jj,jl) 341 ENDIF 342 oa_i(ji,jj,jl) = zindb*zindic*oa_i(ji,jj,jl) 343 344 !--- 2.4 Correction to snow thickness 345 !------------------------------------- 346 zdvres = (zindsn * zindb - 1._wp) * v_s(ji,jj,jl) 347 v_s(ji,jj,jl) = v_s(ji,jj,jl) + zdvres 348 349 !rdm_snw(ji,jj) = rdm_snw(ji,jj) + zdvres * rhosn 350 351 !--- 2.5 Correction to ice thickness 352 !------------------------------------- 353 zdvres = (zindb - 1._wp) * v_i(ji,jj,jl) 354 v_i(ji,jj,jl) = v_i(ji,jj,jl) + zdvres 355 356 !rdm_ice(ji,jj) = rdm_ice(ji,jj) + zdvres * rhoic 357 !sfx_res(ji,jj) = sfx_res(ji,jj) - sm_i(ji,jj,jl) * ( rhoic * zdvres / rdt_ice ) 358 359 !--- 2.6 Snow is transformed into ice if the original ice cover disappears 360 !---------------------------------------------------------------------------- 361 zindg = tms(ji,jj) * MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) ) ) 362 zdvres = zindg * rhosn * v_s(ji,jj,jl) / rau0 363 v_i(ji,jj,jl) = v_i(ji,jj,jl) + zdvres 364 365 zdvres = zindsn*zindb * ( - zindg * v_s(ji,jj,jl) + zindg * v_i(ji,jj,jl) * ( rau0 - rhoic ) / rhosn ) 366 v_s(ji,jj,jl) = v_s(ji,jj,jl) + zdvres 367 368 !--- 2.7 Correction to ice concentrations 369 !-------------------------------------------- 370 a_i(ji,jj,jl) = zindb * a_i(ji,jj,jl) 371 372 !------------------------- 373 ! 2.8) Snow heat content 374 !------------------------- 375 e_s(ji,jj,1,jl) = zindsn * ( MIN ( MAX ( 0.0, e_s(ji,jj,1,jl) ), zbigvalue ) ) 376 377 END DO ! ji 378 END DO ! jj 379 END DO ! jl 380 381 !------------------------ 382 ! 2.9) Ice heat content 383 !------------------------ 384 385 DO jl = 1, jpl 386 DO jk = 1, nlay_i 387 DO jj = 1, jpj 388 DO ji = 1, jpi 389 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) 390 e_i(ji,jj,jk,jl)= zindic * ( MIN ( MAX ( 0.0, e_i(ji,jj,jk,jl) ), zbigvalue ) ) 391 END DO ! ji 392 END DO ! jj 393 END DO !jk 394 END DO !jl 395 396 132 133 !clem debug: it is done in limthd_dh now 134 ! ! Melt of snow 135 ! !-------------- 136 ! DO jl = 1, jpl 137 ! DO jj = 1, jpj 138 ! DO ji = 1, jpi 139 ! IF( v_s(ji,jj,jl) >= epsi20 ) THEN 140 ! ! If snow energy of melting smaller then Lf 141 ! ! Then all snow melts and heat go to the ocean 142 ! !IF ( zEs <= lfus ) THEN 143 ! IF( t_s(ji,jj,1,jl) >= rtt ) THEN 144 ! zdvres = - v_s(ji,jj,jl) 145 ! zEs = - e_s(ji,jj,1,jl) * unit_fac / ( area(ji,jj) * MAX( v_s(ji,jj,jl), epsi20 ) ) ! snow energy of melting (J.m-3) 146 ! ! Contribution to heat flux to the ocean [W.m-2], < 0 147 ! hfx_res(ji,jj) = hfx_res(ji,jj) - zEs * zdvres * r1_rdtice 148 ! ! Contribution to mass flux 149 ! wfx_snw(ji,jj) = wfx_snw(ji,jj) + rhosn * zdvres * r1_rdtice 150 ! ! updates 151 ! v_s (ji,jj,jl) = 0._wp 152 ! ht_s(ji,jj,jl) = 0._wp 153 ! e_s (ji,jj,1,jl) = 0._wp 154 ! t_s (ji,jj,1,jl) = rtt 155 ! ENDIF 156 ! ENDIF 157 ! END DO 158 ! END DO 159 ! END DO 160 !clem debug 161 162 !--- 2.12 Constrain the thickness of the smallest category above 10 cm 163 !---------------------------------------------------------------------- 397 164 DO jm = 1, jpm 398 165 DO jj = 1, jpj 399 166 DO ji = 1, jpi 400 167 jl = ice_cat_bounds(jm,1) 401 !--- 2.12 Constrain the thickness of the smallest category above 5 cm 402 !---------------------------------------------------------------------- 403 zindb = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi10 ) ) 404 ht_i(ji,jj,jl) = zindb*v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl), epsi10) 405 zh = MAX( rone , zindb * hiclim / MAX( ht_i(ji,jj,jl) , epsi10 ) ) 406 ht_s(ji,jj,jl) = ht_s(ji,jj,jl)* zh 407 ht_i(ji,jj,jl) = ht_i(ji,jj,jl)* zh 408 a_i (ji,jj,jl) = a_i(ji,jj,jl) / zh 409 !CLEM 410 v_i (ji,jj,jl) = a_i(ji,jj,jl) * ht_i(ji,jj,jl) 411 v_s (ji,jj,jl) = a_i(ji,jj,jl) * ht_s(ji,jj,jl) 168 IF( v_i(ji,jj,jl) > 0._wp .AND. ht_i(ji,jj,jl) < hiclim ) THEN 169 zh = hiclim / ht_i(ji,jj,jl) 170 ht_s(ji,jj,jl) = ht_s(ji,jj,jl) * zh 171 ht_i(ji,jj,jl) = ht_i(ji,jj,jl) * zh 172 a_i (ji,jj,jl) = a_i(ji,jj,jl) / zh 173 ENDIF 412 174 END DO !ji 413 175 END DO !jj 414 176 END DO !jm 415 177 178 !--- 2.13 ice concentration should not exceed amax 179 !----------------------------------------------------- 416 180 at_i(:,:) = 0.0 417 181 DO jl = 1, jpl 418 182 at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 419 183 END DO 420 421 !--- 2.13 ice concentration should not exceed amax 422 ! (it should not be the case) 423 !----------------------------------------------------- 424 DO jj = 1, jpj 425 DO ji = 1, jpi 426 z_da_ex = MAX( at_i(ji,jj) - amax , 0.0 ) 427 zindb = MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi10 ) ) 428 DO jl = 1, jpl 429 z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi10 ) * zindb 430 a_i(ji,jj,jl) = MAX( 0._wp, a_i(ji,jj,jl) - z_da_i ) 431 ! 432 zinda = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi10 ) ) 433 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zinda 434 !v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! makes ice shrinken but should not be used 184 185 DO jl = 1, jpl 186 DO jj = 1, jpj 187 DO ji = 1, jpi 188 IF( at_i(ji,jj) > amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN 189 a_i(ji,jj,jl) = a_i(ji,jj,jl) * ( 1._wp - ( 1._wp - amax / at_i(ji,jj) ) ) 190 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 191 ENDIF 435 192 END DO 436 193 END DO 437 194 END DO 195 438 196 at_i(:,:) = 0.0 439 197 DO jl = 1, jpl … … 451 209 END DO 452 210 211 ! zap small values 212 !----------------- 213 CALL lim_itd_me_zapsmall 214 453 215 !--------------------- 454 216 ! 2.11) Ice salinity 455 217 !--------------------- 456 ! clem correct bug on smv_i 457 smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:) 458 459 IF ( num_sal == 2 ) THEN ! general case 218 IF ( num_sal == 2 ) THEN 460 219 DO jl = 1, jpl 461 !DO jk = 1, nlay_i462 DO j j = 1, jpj463 DO ji = 1, jpi464 ! salinity stays in bounds465 !clem smv_i(ji,jj,jl) = MAX(MIN((rhoic-rhosn)/rhoic*sss_m(ji,jj),smv_i(ji,jj,jl)),0.1 * v_i(ji,jj,jl) )466 smv_i(ji,jj,jl) = MAX( MIN( s_i_max * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), s_i_min * v_i(ji,jj,jl) )467 i_ice_switch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) ))468 smv_i(ji,jj,jl) = i_ice_switch * smv_i(ji,jj,jl) !+ s_i_min * ( 1._wp - i_ice_switch ) * v_i(ji,jj,jl)469 END DO ! ji470 END DO ! j j471 !END DO !jk220 DO jj = 1, jpj 221 DO ji = 1, jpi 222 zsal = smv_i(ji,jj,jl) 223 smv_i(ji,jj,jl) = sm_i(ji,jj,jl) * v_i(ji,jj,jl) 224 ! salinity stays in bounds 225 i_ice_switch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) ) 226 smv_i(ji,jj,jl) = i_ice_switch * MAX( MIN( s_i_max * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), s_i_min * v_i(ji,jj,jl) ) !+ s_i_min * ( 1._wp - i_ice_switch ) * v_i(ji,jj,jl) 227 ! associated salt flux 228 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice 229 END DO ! ji 230 END DO ! jj 472 231 END DO !jl 473 232 ENDIF 233 474 234 475 235 ! ------------------- … … 501 261 v_ice(:,:) = v_ice(:,:) * tmv(:,:) 502 262 503 !-------------------------------- 504 ! Update mass/salt fluxes (clem) 505 !-------------------------------- 506 DO jl = 1, jpl 507 DO jj = 1, jpj 508 DO ji = 1, jpi 509 diag_res_pr(ji,jj) = diag_res_pr(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) / rdt_ice 510 rdm_ice(ji,jj) = rdm_ice(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * rhoic 511 rdm_snw(ji,jj) = rdm_snw(ji,jj) + ( v_s(ji,jj,jl) - zvsold(ji,jj,jl) ) * rhosn 512 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmvold(ji,jj,jl) ) * rhoic / rdt_ice 513 END DO 514 END DO 515 END DO 263 264 ! ------------------------------------------------- 265 ! Diagnostics 266 ! ------------------------------------------------- 267 d_a_i_thd(:,:,:) = a_i(:,:,:) - old_a_i(:,:,:) 268 d_v_s_thd(:,:,:) = v_s(:,:,:) - old_v_s(:,:,:) 269 d_v_i_thd(:,:,:) = v_i(:,:,:) - old_v_i(:,:,:) 270 d_e_s_thd(:,:,:,:) = e_s(:,:,:,:) - old_e_s(:,:,:,:) 271 d_e_i_thd(:,:,1:nlay_i,:) = e_i(:,:,1:nlay_i,:) - old_e_i(:,:,1:nlay_i,:) 272 !?? d_oa_i_thd(:,:,:) = oa_i (:,:,:) - old_oa_i (:,:,:) 273 d_smv_i_thd(:,:,:) = 0._wp 274 IF( num_sal == 2 ) d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 275 ! diag only (clem) 276 dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday 516 277 517 278 ! ------------------------------- 518 279 !- check conservation (C Rousset) 519 IF (ln_limdiahsb) THEN520 521 zchk_f s = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b522 zchk_f w = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b280 IF( ln_limdiahsb ) THEN 281 zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 282 zchk_fw = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fw_b 283 zchk_ft = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) - zchk_ft_b 523 284 524 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:) , dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice285 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b ) * r1_rdtice - zchk_fw 525 286 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 287 zchk_e_i = glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) * r1_rdtice - zchk_e_i_b * r1_rdtice + zchk_ft 526 288 527 289 zchk_vmin = glob_min(v_i) … … 530 292 531 293 IF(lwp) THEN 532 IF ( ABS( zchk_v_i ) > 1.e- 5 ) WRITE(numout,*) 'violation volume [m3/day] (limupdate2) = ',(zchk_v_i * rday)294 IF ( ABS( zchk_v_i ) > 1.e-4 ) WRITE(numout,*) 'violation volume [kg/day] (limupdate2) = ',(zchk_v_i * rday) 533 295 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limupdate2) = ',(zchk_smv * rday) 296 IF ( ABS( zchk_e_i ) > 1.e-2 ) WRITE(numout,*) 'violation enthalpy [1e9 J] (limupdate2) = ',(zchk_e_i) 534 297 IF ( zchk_vmin < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limupdate2) = ',(zchk_vmin * 1.e-3) 535 298 IF ( zchk_amax > amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax (limupdate2) = ',zchk_amax 536 299 IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limupdate2) = ',zchk_amin 537 300 ENDIF 538 ENDIF301 ENDIF 539 302 !- check conservation (C Rousset) 540 303 ! ------------------------------- … … 596 359 CALL prt_ctl_info(' - Heat / FW fluxes : ') 597 360 CALL prt_ctl_info(' ~~~~~~~~~~~~~~~~~~ ') 598 CALL prt_ctl(tab2d_1=fmmec , clinfo1= ' lim_update2 : fmmec : ', tab2d_2=fhmec , clinfo2= ' fhmec : ')599 361 CALL prt_ctl(tab2d_1=sst_m , clinfo1= ' lim_update2 : sst : ', tab2d_2=sss_m , clinfo2= ' sss : ') 600 CALL prt_ctl(tab2d_1=fhbri , clinfo1= ' lim_update2 : fhbri : ', tab2d_2=fheat_mec , clinfo2= ' fheat_mec : ')601 362 602 363 CALL prt_ctl_info(' ') … … 608 369 ENDIF 609 370 610 CALL wrk_dealloc( jpi,jpj,jpl, internal_melt ) ! integer611 CALL wrk_dealloc( jkmax, zthick0, zqm0 )612 613 CALL wrk_dealloc( jpi,jpj,jpl,zviold, zvsold, zsmvold ) ! clem614 615 371 IF( nn_timing == 1 ) CALL timing_stop('limupdate2') 616 372 END SUBROUTINE lim_update2
Note: See TracChangeset
for help on using the changeset viewer.
