Changeset 4921 for branches/2014/dev_r4650_UKMO13_CICE_changes_take2/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
- Timestamp:
- 2014-11-28T14:59:01+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO13_CICE_changes_take2/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
r4333 r4921 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 ! dummy loop indices 70 INTEGER :: i_ice_switch 71 REAL(wp) :: zh, zsal 72 ! 73 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 92 74 !!------------------------------------------------------------------- 93 75 IF( nn_timing == 1 ) CALL timing_start('limupdate2') 94 76 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 ! ------------------------------- 77 ! conservation test 78 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 79 80 !----------------- 81 ! zap small values 82 !----------------- 83 CALL lim_itd_me_zapsmall 130 84 131 85 CALL lim_var_glo2eqv 132 86 133 !-------------------------------------- 134 ! 2. Review of all pathological cases 135 !-------------------------------------- 136 137 ! clem: useless now 138 !------------------------------------------- 139 ! 2.1) Advection of ice in an ice-free cell 140 !------------------------------------------- 141 ! should be removed since it is treated after dynamics now 142 ! zhimax = 5._wp 143 ! ! first category 144 ! DO jj = 1, jpj 145 ! DO ji = 1, jpi 146 ! !--- the thickness of such an ice is often out of bounds 147 ! !--- thus we recompute a new area while conserving ice volume 148 ! 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 line 153 ! ht_i(ji,jj,1) = hi_max(1) * 0.5_wp 154 ! a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1) 155 ! ENDIF 156 ! END DO 157 ! END DO 158 159 ! zhimax = 20._wp 160 ! ! other categories 161 ! DO jl = 2, jpl 162 ! jm = ice_types(jl) 163 ! DO jj = 1, jpj 164 ! DO ji = 1, jpi 165 ! 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 why 167 ! ! it makes problems when the advected volume and concentration do not seem to be 168 ! ! related with each other 169 ! ! the new thickness is sometimes very big! 170 ! ! and sometimes d_a_i_trp and d_v_i_trp have different sign 171 ! ! which of course is plausible 172 ! ! 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 ) THEN 175 ! 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_wp 176 ! a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl) 177 ! ENDIF 178 ! END DO ! ji 179 ! END DO !jj 180 ! END DO !jl 181 87 !---------------------------------------------------- 88 ! Rebin categories with thickness out of bounds 89 !---------------------------------------------------- 90 IF ( jpl > 1 ) CALL lim_itd_th_reb(1, jpl) 91 92 !---------------------------------------------------------------------- 93 ! Constrain the thickness of the smallest category above hiclim 94 !---------------------------------------------------------------------- 95 DO jj = 1, jpj 96 DO ji = 1, jpi 97 IF( v_i(ji,jj,1) > 0._wp .AND. ht_i(ji,jj,1) < hiclim ) THEN 98 zh = hiclim / ht_i(ji,jj,1) 99 ht_s(ji,jj,1) = ht_s(ji,jj,1) * zh 100 ht_i(ji,jj,1) = ht_i(ji,jj,1) * zh 101 a_i (ji,jj,1) = a_i(ji,jj,1) / zh 102 ENDIF 103 END DO 104 END DO 105 106 !----------------------------------------------------- 107 ! ice concentration should not exceed amax 108 !----------------------------------------------------- 182 109 at_i(:,:) = 0._wp 183 110 DO jl = 1, jpl … … 185 112 END DO 186 113 187 !---------------------------------------------------- 188 ! 2.2) Rebin categories with thickness out of bounds 189 !---------------------------------------------------- 190 DO jm = 1, jpm 191 jbnd1 = ice_cat_bounds(jm,1) 192 jbnd2 = ice_cat_bounds(jm,2) 193 IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_reb(jbnd1, jbnd2, jm) 194 END DO 195 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 114 DO jl = 1, jpl 115 DO jj = 1, jpj 217 116 DO ji = 1, jpi 218 IF( internal_melt(ji,jj,jl) == 1 ) THEN 219 ! initial ice thickness 220 !----------------------- 117 IF( at_i(ji,jj) > amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN 118 a_i(ji,jj,jl) = a_i(ji,jj,jl) * ( 1._wp - ( 1._wp - amax / at_i(ji,jj) ) ) 221 119 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 120 ENDIF 300 121 END DO 301 122 END DO 302 123 END DO 303 304 DO jl = 1, jpl305 DO jj = 1, jpj306 DO ji = 1, jpi307 IF ( internal_melt(ji,jj,jl) == 1 ) THEN308 zdvres = v_s(ji,jj,jl)309 ! release heat310 fheat_res(ji,jj) = fheat_res(ji,jj) + ze_s * zdvres / rdt_ice311 ! release mass312 !rdm_snw(ji,jj) = rdm_snw(ji,jj) - zdvres * rhosn313 !314 v_s(ji,jj,jl) = 0.0315 e_s(ji,jj,1,jl) = 0.0316 ENDIF317 END DO318 END DO319 END DO320 321 zbigvalue = 1.0e+20322 DO jl = 1, jpl323 DO jj = 1, jpj324 DO ji = 1, jpi325 326 !switches327 zindb = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi10 ) )328 !switch = 1 if a_i > 1e-06 and 0 if not329 zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) ) !=1 if hs > 1e-10 and 0 if not330 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) !=1 if hi > 1e-10 and 0 if not331 ! bug fix 25 avril 2007332 zindb = zindb*zindic333 334 !--- 2.3 Correction to ice age335 !------------------------------336 ! IF ((o_i(ji,jj,jl)-1.0)*rday.gt.(rdt_ice*float(numit))) THEN337 ! o_i(ji,jj,jl) = rdt_ice*FLOAT(numit)/rday338 ! ENDIF339 IF ((oa_i(ji,jj,jl)-1.0)*rday.gt.(rdt_ice*numit*a_i(ji,jj,jl))) THEN340 oa_i(ji,jj,jl) = rdt_ice*numit/rday*a_i(ji,jj,jl)341 ENDIF342 oa_i(ji,jj,jl) = zindb*zindic*oa_i(ji,jj,jl)343 344 !--- 2.4 Correction to snow thickness345 !-------------------------------------346 zdvres = (zindsn * zindb - 1._wp) * v_s(ji,jj,jl)347 v_s(ji,jj,jl) = v_s(ji,jj,jl) + zdvres348 349 !rdm_snw(ji,jj) = rdm_snw(ji,jj) + zdvres * rhosn350 351 !--- 2.5 Correction to ice thickness352 !-------------------------------------353 zdvres = (zindb - 1._wp) * v_i(ji,jj,jl)354 v_i(ji,jj,jl) = v_i(ji,jj,jl) + zdvres355 356 !rdm_ice(ji,jj) = rdm_ice(ji,jj) + zdvres * rhoic357 !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 disappears360 !----------------------------------------------------------------------------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) / rau0363 v_i(ji,jj,jl) = v_i(ji,jj,jl) + zdvres364 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) + zdvres367 368 !--- 2.7 Correction to ice concentrations369 !--------------------------------------------370 a_i(ji,jj,jl) = zindb * a_i(ji,jj,jl)371 372 !-------------------------373 ! 2.8) Snow heat content374 !-------------------------375 e_s(ji,jj,1,jl) = zindsn * ( MIN ( MAX ( 0.0, e_s(ji,jj,1,jl) ), zbigvalue ) )376 377 END DO ! ji378 END DO ! jj379 END DO ! jl380 381 !------------------------382 ! 2.9) Ice heat content383 !------------------------384 385 DO jl = 1, jpl386 DO jk = 1, nlay_i387 DO jj = 1, jpj388 DO ji = 1, jpi389 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 ! ji392 END DO ! jj393 END DO !jk394 END DO !jl395 396 397 DO jm = 1, jpm398 DO jj = 1, jpj399 DO ji = 1, jpi400 jl = ice_cat_bounds(jm,1)401 !--- 2.12 Constrain the thickness of the smallest category above 5 cm402 !----------------------------------------------------------------------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)* zh407 ht_i(ji,jj,jl) = ht_i(ji,jj,jl)* zh408 a_i (ji,jj,jl) = a_i(ji,jj,jl) / zh409 !CLEM410 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)412 END DO !ji413 END DO !jj414 END DO !jm415 124 416 125 at_i(:,:) = 0.0 … … 418 127 at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 419 128 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 129 130 ! -------------------------------------- 443 131 ! Final thickness distribution rebinning 444 132 ! -------------------------------------- 445 DO jm = 1, jpm 446 jbnd1 = ice_cat_bounds(jm,1) 447 jbnd2 = ice_cat_bounds(jm,2) 448 IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_reb(jbnd1, jbnd2, jm) 449 IF (ice_ncat_types(jm) .EQ. 1 ) THEN 450 ENDIF 451 END DO 133 IF ( jpl > 1 ) CALL lim_itd_th_reb( 1, jpl ) 134 135 !----------------- 136 ! zap small values 137 !----------------- 138 CALL lim_itd_me_zapsmall 452 139 453 140 !--------------------- 454 141 ! 2.11) Ice salinity 455 142 !--------------------- 456 ! clem correct bug on smv_i 457 smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:) 458 459 IF ( num_sal == 2 ) THEN ! general case 143 IF ( num_sal == 2 ) THEN 460 144 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 !jk145 DO jj = 1, jpj 146 DO ji = 1, jpi 147 zsal = smv_i(ji,jj,jl) 148 smv_i(ji,jj,jl) = sm_i(ji,jj,jl) * v_i(ji,jj,jl) 149 ! salinity stays in bounds 150 i_ice_switch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) ) 151 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) 152 ! associated salt flux 153 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice 154 END DO ! ji 155 END DO ! jj 472 156 END DO !jl 473 157 ENDIF 474 475 ! -------------------476 at_i(:,:) = a_i(:,:,1)477 DO jl = 2, jpl478 at_i(:,:) = a_i(:,:,jl) + at_i(:,:)479 END DO480 158 481 159 !------------------------------------------------------------------------------ … … 486 164 DO jj = 2, jpjm1 487 165 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 side166 IF ( at_i(ji,jj) == 0._wp ) THEN ! what to do if there is no ice 167 IF ( at_i(ji+1,jj) == 0._wp ) u_ice(ji,jj) = 0._wp ! right side 168 IF ( at_i(ji-1,jj) == 0._wp ) u_ice(ji-1,jj) = 0._wp ! left side 169 IF ( at_i(ji,jj+1) == 0._wp ) v_ice(ji,jj) = 0._wp ! upper side 170 IF ( at_i(ji,jj-1) == 0._wp ) v_ice(ji,jj-1) = 0._wp ! bottom side 493 171 ENDIF 494 172 END DO … … 501 179 v_ice(:,:) = v_ice(:,:) * tmv(:,:) 502 180 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 ! ------------------------------- 181 ! ------------------------------------------------- 182 ! Diagnostics 183 ! ------------------------------------------------- 184 d_a_i_thd(:,:,:) = a_i(:,:,:) - a_i_b(:,:,:) 185 d_v_s_thd(:,:,:) = v_s(:,:,:) - v_s_b(:,:,:) 186 d_v_i_thd(:,:,:) = v_i(:,:,:) - v_i_b(:,:,:) 187 d_e_s_thd(:,:,:,:) = e_s(:,:,:,:) - e_s_b(:,:,:,:) 188 d_e_i_thd(:,:,1:nlay_i,:) = e_i(:,:,1:nlay_i,:) - e_i_b(:,:,1:nlay_i,:) 189 !?? d_oa_i_thd(:,:,:) = oa_i (:,:,:) - oa_i_b (:,:,:) 190 d_smv_i_thd(:,:,:) = 0._wp 191 IF( num_sal == 2 ) d_smv_i_thd(:,:,:) = smv_i(:,:,:) - smv_i_b(:,:,:) 192 ! diag only (clem) 193 dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday 194 195 ! heat content variation (W.m-2) 196 DO jj = 1, jpj 197 DO ji = 1, jpi 198 diag_heat_dhc(ji,jj) = ( SUM( d_e_i_trp(ji,jj,1:nlay_i,:) + d_e_i_thd(ji,jj,1:nlay_i,:) ) + & 199 & SUM( d_e_s_trp(ji,jj,1:nlay_s,:) + d_e_s_thd(ji,jj,1:nlay_s,:) ) & 200 & ) * unit_fac * r1_rdtice / area(ji,jj) 201 END DO 202 END DO 203 204 ! conservation test 205 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 541 206 542 207 IF(ln_ctl) THEN ! Control print … … 550 215 CALL prt_ctl(tab2d_1=strength , clinfo1=' lim_update2 : strength :') 551 216 CALL prt_ctl(tab2d_1=u_ice , clinfo1=' lim_update2 : u_ice :', tab2d_2=v_ice , clinfo2=' v_ice :') 552 CALL prt_ctl(tab2d_1= old_u_ice , clinfo1=' lim_update2 : old_u_ice :', tab2d_2=old_v_ice , clinfo2=' old_v_ice:')217 CALL prt_ctl(tab2d_1=u_ice_b , clinfo1=' lim_update2 : u_ice_b :', tab2d_2=v_ice_b , clinfo2=' v_ice_b :') 553 218 554 219 DO jl = 1, jpl … … 563 228 CALL prt_ctl(tab2d_1=o_i (:,:,jl) , clinfo1= ' lim_update2 : o_i : ') 564 229 CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_update2 : a_i : ') 565 CALL prt_ctl(tab2d_1= old_a_i (:,:,jl) , clinfo1= ' lim_update2 : old_a_i: ')230 CALL prt_ctl(tab2d_1=a_i_b (:,:,jl) , clinfo1= ' lim_update2 : a_i_b : ') 566 231 CALL prt_ctl(tab2d_1=d_a_i_thd (:,:,jl) , clinfo1= ' lim_update2 : d_a_i_thd : ') 567 232 CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_update2 : v_i : ') 568 CALL prt_ctl(tab2d_1= old_v_i (:,:,jl) , clinfo1= ' lim_update2 : old_v_i: ')233 CALL prt_ctl(tab2d_1=v_i_b (:,:,jl) , clinfo1= ' lim_update2 : v_i_b : ') 569 234 CALL prt_ctl(tab2d_1=d_v_i_thd (:,:,jl) , clinfo1= ' lim_update2 : d_v_i_thd : ') 570 235 CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_update2 : v_s : ') 571 CALL prt_ctl(tab2d_1= old_v_s (:,:,jl) , clinfo1= ' lim_update2 : old_v_s: ')236 CALL prt_ctl(tab2d_1=v_s_b (:,:,jl) , clinfo1= ' lim_update2 : v_s_b : ') 572 237 CALL prt_ctl(tab2d_1=d_v_s_thd (:,:,jl) , clinfo1= ' lim_update2 : d_v_s_thd : ') 573 238 CALL prt_ctl(tab2d_1=e_i (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2 : e_i1 : ') 574 CALL prt_ctl(tab2d_1= old_e_i (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2 : old_e_i1: ')239 CALL prt_ctl(tab2d_1=e_i_b (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2 : e_i1_b : ') 575 240 CALL prt_ctl(tab2d_1=d_e_i_thd (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2 : de_i1_thd : ') 576 241 CALL prt_ctl(tab2d_1=e_i (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2 : e_i2 : ') 577 CALL prt_ctl(tab2d_1= old_e_i (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2 : old_e_i2: ')242 CALL prt_ctl(tab2d_1=e_i_b (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2 : e_i2_b : ') 578 243 CALL prt_ctl(tab2d_1=d_e_i_thd (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2 : de_i2_thd : ') 579 244 CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_update2 : e_snow : ') 580 CALL prt_ctl(tab2d_1= old_e_s (:,:,1,jl) , clinfo1= ' lim_update2 : old_e_snow: ')245 CALL prt_ctl(tab2d_1=e_s_b (:,:,1,jl) , clinfo1= ' lim_update2 : e_snow_b : ') 581 246 CALL prt_ctl(tab2d_1=d_e_s_thd (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2 : d_e_s_thd : ') 582 247 CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_update2 : smv_i : ') 583 CALL prt_ctl(tab2d_1= old_smv_i (:,:,jl) , clinfo1= ' lim_update2 : old_smv_i: ')248 CALL prt_ctl(tab2d_1=smv_i_b (:,:,jl) , clinfo1= ' lim_update2 : smv_i_b : ') 584 249 CALL prt_ctl(tab2d_1=d_smv_i_thd(:,:,jl) , clinfo1= ' lim_update2 : d_smv_i_thd : ') 585 250 CALL prt_ctl(tab2d_1=oa_i (:,:,jl) , clinfo1= ' lim_update2 : oa_i : ') 586 CALL prt_ctl(tab2d_1=o ld_oa_i (:,:,jl) , clinfo1= ' lim_update2 : old_oa_i: ')251 CALL prt_ctl(tab2d_1=oa_i_b (:,:,jl) , clinfo1= ' lim_update2 : oa_i_b : ') 587 252 CALL prt_ctl(tab2d_1=d_oa_i_thd (:,:,jl) , clinfo1= ' lim_update2 : d_oa_i_thd : ') 588 253 … … 596 261 CALL prt_ctl_info(' - Heat / FW fluxes : ') 597 262 CALL prt_ctl_info(' ~~~~~~~~~~~~~~~~~~ ') 598 CALL prt_ctl(tab2d_1=fmmec , clinfo1= ' lim_update2 : fmmec : ', tab2d_2=fhmec , clinfo2= ' fhmec : ')599 263 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 264 602 265 CALL prt_ctl_info(' ') … … 608 271 ENDIF 609 272 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 273 IF( nn_timing == 1 ) CALL timing_stop('limupdate2') 274 616 275 END SUBROUTINE lim_update2 617 276 #else
Note: See TracChangeset
for help on using the changeset viewer.
