Changeset 4688 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
- Timestamp:
- 2014-06-25T01:39:59+02:00 (10 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90 (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
r4333 r4688 5 5 !!====================================================================== 6 6 !! History : 3.0 ! 2006-04 (M. Vancoppenolle) Original code 7 !! 3.6 ! 2014-06 (C. Rousset) Complete rewriting/cleaning 7 8 !!---------------------------------------------------------------------- 8 9 #if defined key_lim3 … … 39 40 USE lib_fortran ! glob_sum 40 41 USE timing ! Timing 42 USE limcons ! conservation tests 41 43 42 44 IMPLICIT NONE … … 45 47 PUBLIC lim_update2 ! routine called by ice_step 46 48 47 REAL(wp) :: epsi10 = 1.e-10_wp ! - - 48 REAL(wp) :: rzero = 0._wp ! - - 49 REAL(wp) :: rone = 1._wp ! - - 50 49 REAL(wp) :: epsi10 = 1.e-10_wp ! - - 50 REAL(wp) :: epsi20 = 1.e-20_wp 51 51 52 !! * Substitutions 52 53 # include "vectopt_loop_substitute.h90" … … 64 65 !! ** Purpose : Computes update of sea-ice global variables at 65 66 !! the end of the time step. 66 !! Address pathological cases67 !! This place is very important68 !!69 !! ** Method :70 !! Ice speed from ice dynamics71 !! Ice thickness, Snow thickness, Temperatures, Lead fraction72 !! from advection and ice thermodynamics73 67 !! 74 !! ** Action : -75 68 !!--------------------------------------------------------------------- 76 INTEGER :: ji, jj, jk, jl, jm ! dummy loop indices 77 INTEGER :: jbnd1, jbnd2 78 INTEGER :: i_ice_switch 79 INTEGER :: ind_im, layer ! indices for internal melt 80 REAL(wp) :: zweight, zesum, zhimax, z_da_i 81 REAL(wp) :: zinda, zindb, zindsn, zindic 82 REAL(wp) :: zindg, zh, zdvres, zviold2 83 REAL(wp) :: zbigvalue, zvsold2, z_da_ex 84 REAL(wp) :: z_prescr_hi, zat_i_old, ztmelts, ze_s 85 86 INTEGER , POINTER, DIMENSION(:,:,:) :: internal_melt 87 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) 89 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... 69 INTEGER :: ji, jj, jk, jl, jm ! dummy loop indices 70 INTEGER :: jbnd1, jbnd2 71 INTEGER :: i_ice_switch 72 REAL(wp) :: zh, zsal 73 ! 74 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 92 75 !!------------------------------------------------------------------- 93 76 IF( nn_timing == 1 ) CALL timing_start('limupdate2') 94 77 95 CALL wrk_alloc( jpi,jpj,jpl, internal_melt ) ! integer 96 CALL wrk_alloc( jkmax, zthick0, zqm0 ) 97 98 CALL wrk_alloc( jpi,jpj,jpl,zviold, zvsold, zsmvold ) ! clem 99 100 !---------------------------------------------------------------------------------------- 101 ! 1. Computation of trend terms 102 !---------------------------------------------------------------------------------------- 103 !- Trend terms 104 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._wp 111 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 * rday 114 115 ! mass and salt flux init (clem) 116 zviold(:,:,:) = v_i(:,:,:) 117 zvsold(:,:,:) = v_s(:,:,:) 118 zsmvold(:,:,:) = smv_i(:,:,:) 119 120 ! ------------------------------- 121 !- check conservation (C Rousset) 122 IF (ln_limdiahsb) THEN 123 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 124 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(:,:) ) 127 ENDIF 128 !- check conservation (C Rousset) 129 ! ------------------------------- 78 ! conservation test 79 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 80 81 !----------------- 82 ! zap small values 83 !----------------- 84 CALL lim_itd_me_zapsmall 130 85 131 86 CALL lim_var_glo2eqv 132 87 133 !--------------------------------------134 ! 2. Review of all pathological cases135 !--------------------------------------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 at_i(:,:) = 0._wp183 DO jl = 1, jpl184 at_i(:,:) = a_i(:,:,jl) + at_i(:,:)185 END DO186 187 88 !---------------------------------------------------- 188 ! 2.2)Rebin categories with thickness out of bounds89 ! Rebin categories with thickness out of bounds 189 90 !---------------------------------------------------- 190 91 DO jm = 1, jpm … … 194 95 END DO 195 96 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 97 !---------------------------------------------------------------------- 98 ! Constrain the thickness of the smallest category above hiclim 99 !---------------------------------------------------------------------- 397 100 DO jm = 1, jpm 398 101 DO jj = 1, jpj 399 102 DO ji = 1, jpi 400 103 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) 104 IF( v_i(ji,jj,jl) > 0._wp .AND. ht_i(ji,jj,jl) < hiclim ) THEN 105 zh = hiclim / ht_i(ji,jj,jl) 106 ht_s(ji,jj,jl) = ht_s(ji,jj,jl) * zh 107 ht_i(ji,jj,jl) = ht_i(ji,jj,jl) * zh 108 a_i (ji,jj,jl) = a_i(ji,jj,jl) / zh 109 ENDIF 412 110 END DO !ji 413 111 END DO !jj 414 112 END DO !jm 113 114 !----------------------------------------------------- 115 ! ice concentration should not exceed amax 116 !----------------------------------------------------- 117 at_i(:,:) = 0._wp 118 DO jl = 1, jpl 119 at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 120 END DO 121 122 DO jl = 1, jpl 123 DO jj = 1, jpj 124 DO ji = 1, jpi 125 IF( at_i(ji,jj) > amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN 126 a_i(ji,jj,jl) = a_i(ji,jj,jl) * ( 1._wp - ( 1._wp - amax / at_i(ji,jj) ) ) 127 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 128 ENDIF 129 END DO 130 END DO 131 END DO 415 132 416 133 at_i(:,:) = 0.0 … … 418 135 at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 419 136 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 435 END DO 436 END DO 437 END DO 438 at_i(:,:) = 0.0 439 DO jl = 1, jpl 440 at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 441 END DO 442 137 138 ! -------------------------------------- 443 139 ! Final thickness distribution rebinning 444 140 ! -------------------------------------- … … 451 147 END DO 452 148 149 !----------------- 150 ! zap small values 151 !----------------- 152 CALL lim_itd_me_zapsmall 153 453 154 !--------------------- 454 155 ! 2.11) Ice salinity 455 156 !--------------------- 456 ! clem correct bug on smv_i 457 smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:) 458 459 IF ( num_sal == 2 ) THEN ! general case 157 IF ( num_sal == 2 ) THEN 460 158 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 !jk159 DO jj = 1, jpj 160 DO ji = 1, jpi 161 zsal = smv_i(ji,jj,jl) 162 smv_i(ji,jj,jl) = sm_i(ji,jj,jl) * v_i(ji,jj,jl) 163 ! salinity stays in bounds 164 i_ice_switch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) ) 165 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) 166 ! associated salt flux 167 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice 168 END DO ! ji 169 END DO ! jj 472 170 END DO !jl 473 171 ENDIF 474 475 ! -------------------476 at_i(:,:) = a_i(:,:,1)477 DO jl = 2, jpl478 at_i(:,:) = a_i(:,:,jl) + at_i(:,:)479 END DO480 172 481 173 !------------------------------------------------------------------------------ … … 486 178 DO jj = 2, jpjm1 487 179 DO ji = 2, jpim1 488 IF ( at_i(ji,jj) .EQ. 0.0) THEN ! what to do if there is no ice489 IF ( at_i(ji+1,jj) .EQ. 0.0 ) u_ice(ji,jj) = 0.0! right side490 IF ( at_i(ji-1,jj) .EQ. 0.0 ) u_ice(ji-1,jj) = 0.0! left side491 IF ( at_i(ji,jj+1) .EQ. 0.0 ) v_ice(ji,jj) = 0.0! upper side492 IF ( at_i(ji,jj-1) .EQ. 0.0 ) v_ice(ji,jj-1) = 0.0! bottom side180 IF ( at_i(ji,jj) == 0._wp ) THEN ! what to do if there is no ice 181 IF ( at_i(ji+1,jj) == 0._wp ) u_ice(ji,jj) = 0._wp ! right side 182 IF ( at_i(ji-1,jj) == 0._wp ) u_ice(ji-1,jj) = 0._wp ! left side 183 IF ( at_i(ji,jj+1) == 0._wp ) v_ice(ji,jj) = 0._wp ! upper side 184 IF ( at_i(ji,jj-1) == 0._wp ) v_ice(ji,jj-1) = 0._wp ! bottom side 493 185 ENDIF 494 186 END DO … … 501 193 v_ice(:,:) = v_ice(:,:) * tmv(:,:) 502 194 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 516 517 ! ------------------------------- 518 !- check conservation (C Rousset) 519 IF (ln_limdiahsb) THEN 520 521 zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 522 zchk_fw = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 523 524 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice 525 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 526 527 zchk_vmin = glob_min(v_i) 528 zchk_amax = glob_max(SUM(a_i,dim=3)) 529 zchk_amin = glob_min(a_i) 530 531 IF(lwp) THEN 532 IF ( ABS( zchk_v_i ) > 1.e-5 ) WRITE(numout,*) 'violation volume [m3/day] (limupdate2) = ',(zchk_v_i * rday) 533 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limupdate2) = ',(zchk_smv * rday) 534 IF ( zchk_vmin < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limupdate2) = ',(zchk_vmin * 1.e-3) 535 IF ( zchk_amax > amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax (limupdate2) = ',zchk_amax 536 IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limupdate2) = ',zchk_amin 537 ENDIF 538 ENDIF 539 !- check conservation (C Rousset) 540 ! ------------------------------- 195 ! ------------------------------------------------- 196 ! Diagnostics 197 ! ------------------------------------------------- 198 d_a_i_thd(:,:,:) = a_i(:,:,:) - old_a_i(:,:,:) 199 d_v_s_thd(:,:,:) = v_s(:,:,:) - old_v_s(:,:,:) 200 d_v_i_thd(:,:,:) = v_i(:,:,:) - old_v_i(:,:,:) 201 d_e_s_thd(:,:,:,:) = e_s(:,:,:,:) - old_e_s(:,:,:,:) 202 d_e_i_thd(:,:,1:nlay_i,:) = e_i(:,:,1:nlay_i,:) - old_e_i(:,:,1:nlay_i,:) 203 !?? d_oa_i_thd(:,:,:) = oa_i (:,:,:) - old_oa_i (:,:,:) 204 d_smv_i_thd(:,:,:) = 0._wp 205 IF( num_sal == 2 ) d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 206 ! diag only (clem) 207 dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday 208 209 ! heat content variation (W.m-2) 210 DO jj = 1, jpj 211 DO ji = 1, jpi 212 diag_heat_dhc(ji,jj) = ( SUM( d_e_i_trp(ji,jj,1:nlay_i,:) + d_e_i_thd(ji,jj,1:nlay_i,:) ) + & 213 & SUM( d_e_s_trp(ji,jj,1:nlay_s,:) + d_e_s_thd(ji,jj,1:nlay_s,:) ) ) * unit_fac * r1_rdtice / area(ji,jj) 214 END DO 215 END DO 216 217 ! conservation test 218 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 541 219 542 220 IF(ln_ctl) THEN ! Control print … … 596 274 CALL prt_ctl_info(' - Heat / FW fluxes : ') 597 275 CALL prt_ctl_info(' ~~~~~~~~~~~~~~~~~~ ') 598 CALL prt_ctl(tab2d_1=fmmec , clinfo1= ' lim_update2 : fmmec : ', tab2d_2=fhmec , clinfo2= ' fhmec : ')599 276 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 277 602 278 CALL prt_ctl_info(' ') … … 608 284 ENDIF 609 285 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 286 IF( nn_timing == 1 ) CALL timing_stop('limupdate2') 287 616 288 END SUBROUTINE lim_update2 617 289 #else
Note: See TracChangeset
for help on using the changeset viewer.
