- Timestamp:
- 2014-11-27T16:28:53+01:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_CNRS_2014/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90
r4333 r4900 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 … … 32 33 USE par_ice 33 34 USE limitd_th 35 USE limitd_me 34 36 USE limvar 35 37 USE prtctl ! Print control … … 37 39 USE wrk_nemo ! work arrays 38 40 USE lib_fortran ! glob_sum 39 ! Check budget (Rousset)40 41 USE in_out_manager ! I/O manager 41 42 USE iom ! I/O manager 42 43 USE lib_mpp ! MPP library 43 44 USE timing ! Timing 45 USE limcons ! conservation tests 44 46 45 47 IMPLICIT NONE … … 49 51 50 52 REAL(wp) :: epsi10 = 1.e-10_wp ! - - 51 REAL(wp) :: rzero = 0._wp ! - -52 REAL(wp) :: rone = 1._wp ! - -53 53 54 54 !! * Substitutions … … 66 66 !! 67 67 !! ** Purpose : Computes update of sea-ice global variables at 68 !! the end of the time step. 69 !! Address pathological cases 70 !! This place is very important 68 !! the end of the dynamics. 71 69 !! 72 !! ** Method :73 !! Ice speed from ice dynamics74 !! Ice thickness, Snow thickness, Temperatures, Lead fraction75 !! from advection and ice thermodynamics76 !!77 !! ** Action : -78 70 !!--------------------------------------------------------------------- 79 INTEGER :: ji, jj, jk, jl, jm ! dummy loop indices 80 INTEGER :: jbnd1, jbnd2 81 INTEGER :: i_ice_switch 82 INTEGER :: ind_im, layer ! indices for internal melt 83 REAL(wp) :: zweight, zesum, z_da_i, zhimax 84 REAL(wp) :: zinda, zindb, zindsn, zindic 85 REAL(wp) :: zindg, zh, zdvres, zviold2 86 REAL(wp) :: zbigvalue, zvsold2, z_da_ex 87 REAL(wp) :: z_prescr_hi, zat_i_old, ztmelts, ze_s 88 89 REAL(wp), POINTER, DIMENSION(:) :: zthick0, zqm0 ! thickness of the layers and heat contents for 90 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) 91 REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 92 ! mass and salt flux (clem) 93 REAL(wp), POINTER, DIMENSION(:,:,:) :: zviold, zvsold, zsmvold ! old ice volume... 71 INTEGER :: ji, jj, jk, jl, jm ! dummy loop indices 72 INTEGER :: jbnd1, jbnd2 73 INTEGER :: i_ice_switch 74 REAL(wp) :: zsal 75 ! 76 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 94 77 !!------------------------------------------------------------------- 95 78 IF( nn_timing == 1 ) CALL timing_start('limupdate1') 96 79 97 CALL wrk_alloc( jkmax, zthick0, zqm0 ) 98 99 CALL wrk_alloc( jpi,jpj,jpl,zviold, zvsold, zsmvold ) ! clem 100 101 !------------------------------------------------------------------------------ 102 ! 1. Update of Global variables | 103 !------------------------------------------------------------------------------ 104 105 !----------------- 106 ! Trend terms 107 !----------------- 108 d_u_ice_dyn(:,:) = u_ice(:,:) - old_u_ice(:,:) 109 d_v_ice_dyn(:,:) = v_ice(:,:) - old_v_ice(:,:) 110 d_a_i_trp (:,:,:) = a_i (:,:,:) - old_a_i (:,:,:) 111 d_v_s_trp (:,:,:) = v_s (:,:,:) - old_v_s (:,:,:) 112 d_v_i_trp (:,:,:) = v_i (:,:,:) - old_v_i (:,:,:) 113 d_e_s_trp (:,:,:,:) = e_s (:,:,:,:) - old_e_s (:,:,:,:) 114 d_e_i_trp (:,:,:,:) = e_i (:,:,:,:) - old_e_i (:,:,:,:) 115 d_oa_i_trp (:,:,:) = oa_i (:,:,:) - old_oa_i (:,:,:) 116 d_smv_i_trp(:,:,:) = 0._wp 117 IF( num_sal == 2 ) d_smv_i_trp(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 118 119 ! mass and salt flux init (clem) 120 zviold(:,:,:) = v_i(:,:,:) 121 zvsold(:,:,:) = v_s(:,:,:) 122 zsmvold(:,:,:) = smv_i(:,:,:) 123 124 ! ------------------------------- 125 !- check conservation (C Rousset) 126 IF (ln_limdiahsb) THEN 127 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 128 zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 129 zchk_fw_b = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 130 zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 131 ENDIF 132 !- check conservation (C Rousset) 133 ! ------------------------------- 80 IF( ln_limdyn ) THEN 81 82 ! conservation test 83 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate1', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 84 85 !----------------- 86 ! zap small values 87 !----------------- 88 CALL lim_itd_me_zapsmall 134 89 135 90 CALL lim_var_glo2eqv 136 137 !--------------------------------------138 ! 2. Review of all pathological cases139 !--------------------------------------140 141 ! clem: useless now142 !-------------------------------------------143 ! 2.1) Advection of ice in an ice-free cell144 !-------------------------------------------145 ! should be removed since it is treated after dynamics now146 ! zhimax = 5._wp147 ! ! first category148 ! DO jj = 1, jpj149 ! DO ji = 1, jpi150 ! !--- the thickness of such an ice is often out of bounds151 ! !--- thus we recompute a new area while conserving ice volume152 ! zat_i_old = SUM( old_a_i(ji,jj,:) )153 ! zindb = MAX( 0._wp, SIGN( 1._wp, ABS( d_a_i_trp(ji,jj,1) ) - epsi10 ) )154 ! IF( ( ABS( d_v_i_trp(ji,jj,1) ) / MAX( ABS( d_a_i_trp(ji,jj,1) ), epsi10 ) * zindb .GT. zhimax ) &155 ! & .AND.( ( v_i(ji,jj,1) / MAX( a_i(ji,jj,1), epsi10 ) * zindb ) .GT. zhimax ) &156 ! & .AND.( zat_i_old .LT. 1.e-6 ) ) THEN ! new line157 ! ht_i(ji,jj,1) = hi_max(1) * 0.5_wp158 ! a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1)159 ! ENDIF160 ! END DO161 ! END DO162 !163 ! zhimax = 20._wp164 ! ! other categories165 ! DO jl = 2, jpl166 ! jm = ice_types(jl)167 ! DO jj = 1, jpj168 ! DO ji = 1, jpi169 ! zindb = MAX( rzero, SIGN( rone, ABS( d_a_i_trp(ji,jj,jl) ) - epsi10 ) )170 ! ! this correction is very tricky... sometimes, advection gets wrong i don't know why171 ! ! it makes problems when the advected volume and concentration do not seem to be172 ! ! related with each other173 ! ! the new thickness is sometimes very big!174 ! ! and sometimes d_a_i_trp and d_v_i_trp have different sign175 ! ! which of course is plausible176 ! ! but fuck! it fucks everything up :)177 ! IF ( ( ABS( d_v_i_trp(ji,jj,jl) ) / MAX( ABS( d_a_i_trp(ji,jj,jl) ), epsi10 ) * zindb .GT. zhimax ) &178 ! & .AND. ( v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zindb ) .GT. zhimax ) THEN179 ! 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_wp180 ! a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl)181 ! ENDIF182 ! END DO ! ji183 ! END DO !jj184 ! END DO !jl185 91 186 at_i(:,:) = 0._wp187 DO jl = 1, jpl188 at_i(:,:) = a_i(:,:,jl) + at_i(:,:)189 END DO190 191 92 !---------------------------------------------------- 192 ! 2.2)Rebin categories with thickness out of bounds93 ! Rebin categories with thickness out of bounds 193 94 !---------------------------------------------------- 194 95 DO jm = 1, jpm … … 203 104 END DO 204 105 205 zbigvalue = 1.0e+20 206 207 DO jl = 1, jpl 208 DO jj = 1, jpj 106 !---------------------------------------------------- 107 ! ice concentration should not exceed amax 108 !----------------------------------------------------- 109 DO jl = 1, jpl 110 DO jj = 1, jpj 209 111 DO ji = 1, jpi 210 211 !switches 212 zindb = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi10 ) ) 213 !switch = 1 if a_i > 1e-06 and 0 if not 214 zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) ) !=1 if hs > 1e-10 and 0 if not 215 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) !=1 if hi > 1e-10 and 0 if not 216 ! bug fix 25 avril 2007 217 zindb = zindb*zindic 218 219 !--- 2.3 Correction to ice age 220 !------------------------------ 221 ! IF ((o_i(ji,jj,jl)-1.0)*rday.gt.(rdt_ice*float(numit))) THEN 222 ! o_i(ji,jj,jl) = rdt_ice*FLOAT(numit)/rday 223 ! ENDIF 224 IF ((oa_i(ji,jj,jl)-1.0)*rday.gt.(rdt_ice*numit*a_i(ji,jj,jl))) THEN 225 oa_i(ji,jj,jl) = rdt_ice*numit/rday*a_i(ji,jj,jl) 112 IF( at_i(ji,jj) > amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN 113 a_i(ji,jj,jl) = a_i(ji,jj,jl) * ( 1._wp - ( 1._wp - amax / at_i(ji,jj) ) ) 114 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 226 115 ENDIF 227 oa_i(ji,jj,jl) = zindb*zindic*oa_i(ji,jj,jl) 228 229 !--- 2.4 Correction to snow thickness 230 !------------------------------------- 231 ! ! snow thickness has to be greater than 0, and if ice concentration smaller than 1e-6 then hs = 0 232 ! v_s(ji,jj,jl) = MAX( zindb * v_s(ji,jj,jl), 0.0) 233 ! snow thickness cannot be smaller than 1e-6 234 zdvres = (zindsn * zindb - 1._wp) * v_s(ji,jj,jl) 235 v_s(ji,jj,jl) = v_s(ji,jj,jl) + zdvres 236 237 !rdm_snw(ji,jj) = rdm_snw(ji,jj) + zdvres * rhosn 238 239 !--- 2.5 Correction to ice thickness 240 !------------------------------------- 241 zdvres = (zindb - 1._wp) * v_i(ji,jj,jl) 242 v_i(ji,jj,jl) = v_i(ji,jj,jl) + zdvres 243 244 !rdm_ice(ji,jj) = rdm_ice(ji,jj) + zdvres * rhoic 245 !sfx_res(ji,jj) = sfx_res(ji,jj) - sm_i(ji,jj,jl) * ( rhoic * zdvres / rdt_ice ) 246 247 !--- 2.6 Snow is transformed into ice if the original ice cover disappears 248 !---------------------------------------------------------------------------- 249 zindg = tms(ji,jj) * MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) ) ) 250 zdvres = zindg * rhosn * v_s(ji,jj,jl) / rau0 251 v_i(ji,jj,jl) = v_i(ji,jj,jl) + zdvres 252 253 zdvres = zindsn*zindb * ( - zindg * v_s(ji,jj,jl) + zindg * v_i(ji,jj,jl) * ( rau0 - rhoic ) / rhosn ) 254 v_s(ji,jj,jl) = v_s(ji,jj,jl) + zdvres 255 256 !--- 2.7 Correction to ice concentrations 257 !-------------------------------------------- 258 ! if greater than 0, ice concentration cannot be smaller than 1e-10 259 a_i(ji,jj,jl) = zindb * a_i(ji,jj,jl) 260 261 !------------------------- 262 ! 2.8) Snow heat content 263 !------------------------- 264 e_s(ji,jj,1,jl) = zindsn * ( MIN ( MAX ( 0._wp, e_s(ji,jj,1,jl) ), zbigvalue ) ) 265 266 END DO ! ji 267 END DO ! jj 268 END DO ! jl 269 270 !------------------------ 271 ! 2.9) Ice heat content 272 !------------------------ 273 274 DO jl = 1, jpl 275 DO jk = 1, nlay_i 276 DO jj = 1, jpj 277 DO ji = 1, jpi 278 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) 279 e_i(ji,jj,jk,jl)= zindic * ( MIN ( MAX ( 0.0, e_i(ji,jj,jk,jl) ), zbigvalue ) ) 280 END DO ! ji 281 END DO ! jj 282 END DO !jk 283 END DO !jl 284 116 END DO 117 END DO 118 END DO 119 285 120 at_i(:,:) = 0._wp 286 121 DO jl = 1, jpl 287 122 at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 288 123 END DO 289 290 !--- 2.13 ice concentration should not exceed amax 291 ! (it should not be the case) 292 !----------------------------------------------------- 293 DO jj = 1, jpj 294 DO ji = 1, jpi 295 z_da_ex = MAX( at_i(ji,jj) - amax , 0.0 ) 296 zindb = MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi10 ) ) 297 DO jl = 1, jpl 298 z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi10 ) * zindb 299 a_i(ji,jj,jl) = MAX( 0._wp, a_i(ji,jj,jl) - z_da_i ) 300 ! 301 zinda = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi10 ) ) 302 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zinda 303 !v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! makes ice shrinken but should not be used 304 END DO 305 END DO 306 END DO 307 at_i(:,:) = a_i(:,:,1) 308 DO jl = 2, jpl 309 at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 310 END DO 311 312 124 125 ! -------------------------------------- 313 126 ! Final thickness distribution rebinning 314 127 ! -------------------------------------- … … 321 134 END DO 322 135 136 !----------------- 137 ! zap small values 138 !----------------- 139 CALL lim_itd_me_zapsmall 323 140 324 141 !--------------------- 325 ! 2.11) Ice salinity142 ! Ice salinity bounds 326 143 !--------------------- 327 ! clem correct bug on smv_i 328 smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:) 329 330 IF ( num_sal == 2 ) THEN ! general case 144 IF ( num_sal == 2 ) THEN 331 145 DO jl = 1, jpl 332 !DO jk = 1, nlay_i 333 DO jj = 1, jpj 334 DO ji = 1, jpi 335 ! salinity stays in bounds 336 !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) ) 337 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) ) 338 i_ice_switch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) ) ) 339 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) 340 END DO ! ji 341 END DO ! jj 342 !END DO !jk 343 END DO !jl 344 ENDIF 345 346 at_i(:,:) = a_i(:,:,1) 347 DO jl = 2, jpl 348 at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 349 END DO 350 351 352 !-------------------------------- 353 ! Update mass/salt fluxes (clem) 354 !-------------------------------- 355 DO jl = 1, jpl 356 DO jj = 1, jpj 357 DO ji = 1, jpi 358 diag_res_pr(ji,jj) = diag_res_pr(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) / rdt_ice 359 rdm_ice(ji,jj) = rdm_ice(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * rhoic 360 rdm_snw(ji,jj) = rdm_snw(ji,jj) + ( v_s(ji,jj,jl) - zvsold(ji,jj,jl) ) * rhosn 361 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmvold(ji,jj,jl) ) * rhoic / rdt_ice 146 DO jj = 1, jpj 147 DO ji = 1, jpi 148 zsal = smv_i(ji,jj,jl) 149 smv_i(ji,jj,jl) = sm_i(ji,jj,jl) * v_i(ji,jj,jl) 150 ! salinity stays in bounds 151 i_ice_switch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) ) 152 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) ) 153 ! associated salt flux 154 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice 155 END DO 362 156 END DO 363 157 END DO 364 END DO365 366 ! -------------------------------367 !- check conservation (C Rousset)368 IF (ln_limdiahsb) THEN369 370 zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b371 zchk_fw = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b372 373 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice374 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic )375 376 zchk_vmin = glob_min(v_i)377 zchk_amax = glob_max(SUM(a_i,dim=3))378 zchk_amin = glob_min(a_i)379 380 IF(lwp) THEN381 IF ( ABS( zchk_v_i ) > 1.e-5 ) WRITE(numout,*) 'violation volume [m3/day] (limupdate1) = ',(zchk_v_i * rday)382 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limupdate1) = ',(zchk_smv * rday)383 IF ( zchk_vmin < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limupdate1) = ',(zchk_vmin * 1.e-3)384 IF ( zchk_amax > amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax (limupdate1) = ',zchk_amax385 IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limupdate1) = ',zchk_amin386 ENDIF387 158 ENDIF 388 !- check conservation (C Rousset) 389 ! ------------------------------- 159 160 ! ------------------------------------------------- 161 ! Diagnostics 162 ! ------------------------------------------------- 163 d_u_ice_dyn(:,:) = u_ice(:,:) - old_u_ice(:,:) 164 d_v_ice_dyn(:,:) = v_ice(:,:) - old_v_ice(:,:) 165 d_a_i_trp (:,:,:) = a_i (:,:,:) - old_a_i (:,:,:) 166 d_v_s_trp (:,:,:) = v_s (:,:,:) - old_v_s (:,:,:) 167 d_v_i_trp (:,:,:) = v_i (:,:,:) - old_v_i (:,:,:) 168 d_e_s_trp (:,:,:,:) = e_s (:,:,:,:) - old_e_s (:,:,:,:) 169 d_e_i_trp (:,:,1:nlay_i,:) = e_i (:,:,1:nlay_i,:) - old_e_i(:,:,1:nlay_i,:) 170 d_oa_i_trp (:,:,:) = oa_i (:,:,:) - old_oa_i (:,:,:) 171 d_smv_i_trp(:,:,:) = 0._wp 172 IF( num_sal == 2 ) d_smv_i_trp(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 173 174 ! conservation test 175 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limupdate1', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 390 176 391 177 IF(ln_ctl) THEN ! Control print … … 446 232 CALL prt_ctl_info(' - Heat / FW fluxes : ') 447 233 CALL prt_ctl_info(' ~~~~~~~~~~~~~~~~~~ ') 448 CALL prt_ctl(tab2d_1=fmmec , clinfo1= ' lim_update1 : fmmec : ', tab2d_2=fhmec , clinfo2= ' fhmec : ')449 234 CALL prt_ctl(tab2d_1=sst_m , clinfo1= ' lim_update1 : sst : ', tab2d_2=sss_m , clinfo2= ' sss : ') 450 CALL prt_ctl(tab2d_1=fhbri , clinfo1= ' lim_update1 : fhbri : ', tab2d_2=fheat_mec , clinfo2= ' fheat_mec : ')451 235 452 236 CALL prt_ctl_info(' ') … … 458 242 ENDIF 459 243 460 461 CALL wrk_dealloc( jkmax, zthick0, zqm0 ) 462 463 CALL wrk_dealloc( jpi,jpj,jpl,zviold, zvsold, zsmvold ) ! clem 244 ENDIF ! ln_limdyn 464 245 465 246 IF( nn_timing == 1 ) CALL timing_stop('limupdate1')
Note: See TracChangeset
for help on using the changeset viewer.
