Changeset 5260 for branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.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/limupdate2.F90
- Property svn:keywords set to Id
r4333 r5260 5 5 !!====================================================================== 6 6 !! History : 3.0 ! 2006-04 (M. Vancoppenolle) Original code 7 !! 3.5 ! 2014-06 (C. Rousset) Complete rewriting/cleaning 7 8 !!---------------------------------------------------------------------- 8 9 #if defined key_lim3 … … 12 13 !! lim_update2 : computes update of sea-ice global variables from trend terms 13 14 !!---------------------------------------------------------------------- 14 USE limrhg ! ice rheology15 16 USE dom_oce17 USE oce ! dynamics and tracers variables18 USE in_out_manager19 15 USE sbc_oce ! Surface boundary condition: ocean fields 20 16 USE sbc_ice ! Surface boundary condition: ice fields 21 17 USE dom_ice 18 USE dom_oce 22 19 USE phycst ! physical constants 23 20 USE ice 24 USE limdyn25 USE limtrp26 USE limthd27 USE limsbc28 USE limdiahsb29 USE limwri30 USE limrst31 21 USE thd_ice ! LIM thermodynamic sea-ice variables 32 USE par_ice33 22 USE limitd_th 34 USE limitd_me35 23 USE limvar 36 USE prtctl ! Print control 37 USE lbclnk ! lateral boundary condition - MPP exchanges 38 USE wrk_nemo ! work arrays 39 USE lib_fortran ! glob_sum 24 USE prtctl ! Print control 25 USE lbclnk ! lateral boundary condition - MPP exchanges 26 USE wrk_nemo ! work arrays 40 27 USE timing ! Timing 28 USE limcons ! conservation tests 29 USE limctl 30 USE lib_mpp ! MPP library 31 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 32 USE in_out_manager 41 33 42 34 IMPLICIT NONE … … 45 37 PUBLIC lim_update2 ! routine called by ice_step 46 38 47 REAL(wp) :: epsi10 = 1.e-10_wp ! - -48 REAL(wp) :: rzero = 0._wp ! - -49 REAL(wp) :: rone = 1._wp ! - -50 51 39 !! * Substitutions 52 40 # include "vectopt_loop_substitute.h90" 53 41 !!---------------------------------------------------------------------- 54 42 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 55 !! $Id : limupdate.F90 3294 2012-01-28 16:44:18Z rblod$43 !! $Id$ 56 44 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 57 45 !!---------------------------------------------------------------------- 58 46 CONTAINS 59 47 60 SUBROUTINE lim_update2 48 SUBROUTINE lim_update2( kt ) 61 49 !!------------------------------------------------------------------- 62 50 !! *** ROUTINE lim_update2 *** … … 64 52 !! ** Purpose : Computes update of sea-ice global variables at 65 53 !! 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 54 !! 74 !! ** Action : -75 55 !!--------------------------------------------------------------------- 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... 56 INTEGER, INTENT(in) :: kt ! number of iteration 57 INTEGER :: ji, jj, jk, jl ! dummy loop indices 58 REAL(wp) :: zsal 59 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 92 60 !!------------------------------------------------------------------- 93 61 IF( nn_timing == 1 ) CALL timing_start('limupdate2') 94 62 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(:,:) ) 63 IF( kt == nit000 .AND. lwp ) THEN 64 WRITE(numout,*) ' lim_update2 ' 65 WRITE(numout,*) ' ~~~~~~~~~~~ ' 127 66 ENDIF 128 !- check conservation (C Rousset) 129 ! ------------------------------- 130 131 CALL lim_var_glo2eqv 132 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 67 68 ! conservation test 69 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 70 71 !---------------------------------------------------------------------- 72 ! Constrain the thickness of the smallest category above himin 73 !---------------------------------------------------------------------- 74 DO jj = 1, jpj 75 DO ji = 1, jpi 76 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,1) - epsi20 ) ) !0 if no ice and 1 if yes 77 ht_i(ji,jj,1) = v_i (ji,jj,1) / MAX( a_i(ji,jj,1) , epsi20 ) * rswitch 78 IF( v_i(ji,jj,1) > 0._wp .AND. ht_i(ji,jj,1) < rn_himin ) THEN 79 a_i (ji,jj,1) = a_i (ji,jj,1) * ht_i(ji,jj,1) / rn_himin 80 oa_i(ji,jj,1) = oa_i(ji,jj,1) * ht_i(ji,jj,1) / rn_himin 81 ENDIF 82 END DO 83 END DO 84 85 !----------------------------------------------------- 86 ! ice concentration should not exceed amax 87 !----------------------------------------------------- 182 88 at_i(:,:) = 0._wp 183 89 DO jl = 1, jpl … … 185 91 END DO 186 92 187 !----------------------------------------------------188 ! 2.2) Rebin categories with thickness out of bounds189 !----------------------------------------------------190 DO jm = 1, jpm191 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 DO195 196 !---------------------------------197 ! 2.3) Melt of an internal layer 198 !--------------------- ------------199 internal_melt(:,:,:) = 0200 201 DO jl = 1, jpl202 DO j k = 1, nlay_i93 DO jl = 1, jpl 94 DO jj = 1, jpj 95 DO ji = 1, jpi 96 IF( at_i(ji,jj) > rn_amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN 97 a_i (ji,jj,jl) = a_i (ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) ) 98 oa_i(ji,jj,jl) = oa_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) ) 99 ENDIF 100 END DO 101 END DO 102 END DO 103 104 !--------------------- 105 ! Ice salinity 106 !--------------------- 107 IF ( nn_icesal == 2 ) THEN 108 DO jl = 1, jpl 203 109 DO jj = 1, jpj 204 110 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 111 zsal = smv_i(ji,jj,jl) 112 ! salinity stays in bounds 113 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) ) 114 smv_i(ji,jj,jl) = rswitch * MAX( MIN( rn_simax * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), rn_simin * v_i(ji,jj,jl) ) 115 ! associated salt flux 116 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice 117 END DO 300 118 END DO 301 119 END DO 302 END DO303 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 416 at_i(:,:) = 0.0417 DO jl = 1, jpl418 at_i(:,:) = a_i(:,:,jl) + at_i(:,:)419 END DO420 421 !--- 2.13 ice concentration should not exceed amax422 ! (it should not be the case)423 !-----------------------------------------------------424 DO jj = 1, jpj425 DO ji = 1, jpi426 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, jpl429 z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi10 ) * zindb430 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 ) * zinda434 !v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! makes ice shrinken but should not be used435 END DO436 END DO437 END DO438 at_i(:,:) = 0.0439 DO jl = 1, jpl440 at_i(:,:) = a_i(:,:,jl) + at_i(:,:)441 END DO442 443 ! Final thickness distribution rebinning444 ! --------------------------------------445 DO jm = 1, jpm446 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 ) THEN450 ENDIF451 END DO452 453 !---------------------454 ! 2.11) Ice salinity455 !---------------------456 ! clem correct bug on smv_i457 smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:)458 459 IF ( num_sal == 2 ) THEN ! general case460 DO jl = 1, jpl461 !DO jk = 1, nlay_i462 DO jj = 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 ! jj471 !END DO !jk472 END DO !jl473 120 ENDIF 474 121 475 ! ------------------- 476 at_i(:,:) = a_i(:,:,1) 477 DO jl = 2, jpl 478 at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 479 END DO 122 !---------------------------------------------------- 123 ! Rebin categories with thickness out of bounds 124 !---------------------------------------------------- 125 IF ( jpl > 1 ) CALL lim_itd_th_reb( 1, jpl ) 126 127 !----------------- 128 ! zap small values 129 !----------------- 130 CALL lim_var_zapsmall 480 131 481 132 !------------------------------------------------------------------------------ 482 ! 2)Corrections to avoid wrong values |133 ! Corrections to avoid wrong values | 483 134 !------------------------------------------------------------------------------ 484 135 ! Ice drift … … 486 137 DO jj = 2, jpjm1 487 138 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 side139 IF ( at_i(ji,jj) == 0._wp ) THEN ! what to do if there is no ice 140 IF ( at_i(ji+1,jj) == 0._wp ) u_ice(ji,jj) = 0._wp ! right side 141 IF ( at_i(ji-1,jj) == 0._wp ) u_ice(ji-1,jj) = 0._wp ! left side 142 IF ( at_i(ji,jj+1) == 0._wp ) v_ice(ji,jj) = 0._wp ! upper side 143 IF ( at_i(ji,jj-1) == 0._wp ) v_ice(ji,jj-1) = 0._wp ! bottom side 493 144 ENDIF 494 145 END DO … … 498 149 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 499 150 !mask velocities 500 u_ice(:,:) = u_ice(:,:) * tmu(:,:)501 v_ice(:,:) = v_ice(:,:) * tmv(:,:)151 u_ice(:,:) = u_ice(:,:) * umask(:,:,1) 152 v_ice(:,:) = v_ice(:,:) * vmask(:,:,1) 502 153 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 ! ------------------------------- 154 ! ------------------------------------------------- 155 ! Diagnostics 156 ! ------------------------------------------------- 157 DO jl = 1, jpl 158 oa_i(:,:,jl) = oa_i(:,:,jl) + a_i(:,:,jl) * rdt_ice / rday ! ice natural aging 159 afx_thd(:,:) = afx_thd(:,:) + ( a_i(:,:,jl) - a_i_b(:,:,jl) ) * r1_rdtice 160 END DO 161 afx_tot = afx_thd + afx_dyn 162 163 DO jj = 1, jpj 164 DO ji = 1, jpi 165 ! heat content variation (W.m-2) 166 diag_heat(ji,jj) = diag_heat(ji,jj) - & 167 & ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) + & 168 & SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) ) & 169 & ) * r1_rdtice 170 ! salt, volume 171 diag_smvi(ji,jj) = diag_smvi(ji,jj) + SUM( smv_i(ji,jj,:) - smv_i_b(ji,jj,:) ) * rhoic * r1_rdtice 172 diag_vice(ji,jj) = diag_vice(ji,jj) + SUM( v_i (ji,jj,:) - v_i_b (ji,jj,:) ) * rhoic * r1_rdtice 173 diag_vsnw(ji,jj) = diag_vsnw(ji,jj) + SUM( v_s (ji,jj,:) - v_s_b (ji,jj,:) ) * rhosn * r1_rdtice 174 END DO 175 END DO 176 177 ! conservation test 178 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 179 180 ! for outputs 181 CALL lim_var_glo2eqv 182 CALL lim_var_agg(2) 183 184 ! ------------------------------------------------- 185 ! control prints 186 ! ------------------------------------------------- 187 IF( ln_icectl ) CALL lim_prt( kt, iiceprt, jiceprt, 2, ' - Final state - ' ) ! control print 541 188 542 189 IF(ln_ctl) THEN ! Control print … … 544 191 CALL prt_ctl_info(' - Cell values : ') 545 192 CALL prt_ctl_info(' ~~~~~~~~~~~~~ ') 546 CALL prt_ctl(tab2d_1= area, clinfo1=' lim_update2 : cell area :')193 CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_update2 : cell area :') 547 194 CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_update2 : at_i :') 548 195 CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_update2 : vt_i :') … … 550 197 CALL prt_ctl(tab2d_1=strength , clinfo1=' lim_update2 : strength :') 551 198 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:')199 CALL prt_ctl(tab2d_1=u_ice_b , clinfo1=' lim_update2 : u_ice_b :', tab2d_2=v_ice_b , clinfo2=' v_ice_b :') 553 200 554 201 DO jl = 1, jpl … … 563 210 CALL prt_ctl(tab2d_1=o_i (:,:,jl) , clinfo1= ' lim_update2 : o_i : ') 564 211 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 : ') 566 CALL prt_ctl(tab2d_1=d_a_i_thd (:,:,jl) , clinfo1= ' lim_update2 : d_a_i_thd : ') 212 CALL prt_ctl(tab2d_1=a_i_b (:,:,jl) , clinfo1= ' lim_update2 : a_i_b : ') 567 213 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 : ') 569 CALL prt_ctl(tab2d_1=d_v_i_thd (:,:,jl) , clinfo1= ' lim_update2 : d_v_i_thd : ') 214 CALL prt_ctl(tab2d_1=v_i_b (:,:,jl) , clinfo1= ' lim_update2 : v_i_b : ') 570 215 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 : ') 572 CALL prt_ctl(tab2d_1=d_v_s_thd (:,:,jl) , clinfo1= ' lim_update2 : d_v_s_thd : ') 573 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 : ') 575 CALL prt_ctl(tab2d_1=d_e_i_thd (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2 : de_i1_thd : ') 576 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 : ') 578 CALL prt_ctl(tab2d_1=d_e_i_thd (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2 : de_i2_thd : ') 216 CALL prt_ctl(tab2d_1=v_s_b (:,:,jl) , clinfo1= ' lim_update2 : v_s_b : ') 217 CALL prt_ctl(tab2d_1=e_i (:,:,1,jl) , clinfo1= ' lim_update2 : e_i1 : ') 218 CALL prt_ctl(tab2d_1=e_i_b (:,:,1,jl) , clinfo1= ' lim_update2 : e_i1_b : ') 219 CALL prt_ctl(tab2d_1=e_i (:,:,2,jl) , clinfo1= ' lim_update2 : e_i2 : ') 220 CALL prt_ctl(tab2d_1=e_i_b (:,:,2,jl) , clinfo1= ' lim_update2 : e_i2_b : ') 579 221 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 : ') 581 CALL prt_ctl(tab2d_1=d_e_s_thd (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2 : d_e_s_thd : ') 222 CALL prt_ctl(tab2d_1=e_s_b (:,:,1,jl) , clinfo1= ' lim_update2 : e_snow_b : ') 582 223 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 : ') 584 CALL prt_ctl(tab2d_1=d_smv_i_thd(:,:,jl) , clinfo1= ' lim_update2 : d_smv_i_thd : ') 224 CALL prt_ctl(tab2d_1=smv_i_b (:,:,jl) , clinfo1= ' lim_update2 : smv_i_b : ') 585 225 CALL prt_ctl(tab2d_1=oa_i (:,:,jl) , clinfo1= ' lim_update2 : oa_i : ') 586 CALL prt_ctl(tab2d_1=old_oa_i (:,:,jl) , clinfo1= ' lim_update2 : old_oa_i : ') 587 CALL prt_ctl(tab2d_1=d_oa_i_thd (:,:,jl) , clinfo1= ' lim_update2 : d_oa_i_thd : ') 226 CALL prt_ctl(tab2d_1=oa_i_b (:,:,jl) , clinfo1= ' lim_update2 : oa_i_b : ') 588 227 589 228 DO jk = 1, nlay_i … … 596 235 CALL prt_ctl_info(' - Heat / FW fluxes : ') 597 236 CALL prt_ctl_info(' ~~~~~~~~~~~~~~~~~~ ') 598 CALL prt_ctl(tab2d_1=fmmec , clinfo1= ' lim_update2 : fmmec : ', tab2d_2=fhmec , clinfo2= ' fhmec : ')599 237 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 238 602 239 CALL prt_ctl_info(' ') … … 608 245 ENDIF 609 246 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 247 IF( nn_timing == 1 ) CALL timing_stop('limupdate2') 248 616 249 END SUBROUTINE lim_update2 617 250 #else
Note: See TracChangeset
for help on using the changeset viewer.