- Timestamp:
- 2013-07-09T17:41:20+02:00 (11 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_r3406_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90
r3938 r3963 86 86 INTEGER :: i_ice_switch 87 87 INTEGER :: ind_im, layer ! indices for internal melt 88 REAL(wp) :: zweight, zesum, z_da_i 89 REAL(wp) :: zind b, zindsn, zindic88 REAL(wp) :: zweight, zesum, z_da_i, zhimax 89 REAL(wp) :: zinda, zindb, zindsn, zindic 90 90 REAL(wp) :: zindg, zh, zdvres, zviold2 91 91 REAL(wp) :: zbigvalue, zvsold2, z_da_ex … … 94 94 REAL(wp), POINTER, DIMENSION(:) :: zthick0, zqm0 ! thickness of the layers and heat contents for 95 95 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) 96 REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 96 97 ! mass and salt flux (clem) 97 98 REAL(wp), POINTER, DIMENSION(:,:,:) :: zviold, zvsold, zsmvold ! old ice volume... … … 182 183 ! 2. Review of all pathological cases 183 184 !-------------------------------------- 185 186 !------------------------------------------- 187 ! 2.1) Advection of ice in an ice-free cell 188 !------------------------------------------- 189 ! should be removed since it is treated after dynamics now 190 191 ! !IF( ln_nicep ) THEN 192 ! WRITE(numout,*) ' limupdate1 - before h correction ' 193 ! WRITE(numout,*) ' a_i : ', a_i(12, 44, 1:jpl) 194 ! WRITE(numout,*) ' at_i : ', at_i(12,44) 195 ! WRITE(numout,*) ' v_i : ', v_i(12, 44, 1:jpl) 196 ! WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl) 197 ! !ENDIF 198 ! 199 zhimax = 1._wp 200 ! first category 201 DO jj = 1, jpj 202 DO ji = 1, jpi 203 !--- the thickness of such an ice is often out of bounds 204 !--- thus we recompute a new area while conserving ice volume 205 zat_i_old = SUM( old_a_i(ji,jj,:) ) 206 zindb = MAX( 0._wp, SIGN( 1._wp, ABS( d_a_i_trp(ji,jj,1)) - epsi10 ) ) 207 IF ( ( ABS(d_v_i_trp(ji,jj,1))/MAX(ABS(d_a_i_trp(ji,jj,1)),epsi10)*zindb .GT. zhimax) & 208 .AND.( ( v_i(ji,jj,1)/MAX(a_i(ji,jj,1),epsi10)*zindb) .GT. zhimax ) & 209 .AND.( zat_i_old .LT. 1.e-6 ) ) THEN ! new line 210 ht_i(ji,jj,1) = hi_max(1) / 2.0 211 a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1) 212 ENDIF 213 END DO 214 END DO 215 216 217 ! !IF( ln_nicep ) THEN 218 ! at_i(:,:) = 0._wp 219 ! DO jl = 1, jpl 220 ! at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 221 ! END DO 222 ! WRITE(numout,*) ' limupdate1 - after h correction 1 ' 223 ! WRITE(numout,*) ' a_i : ', a_i(12, 44, 1:jpl) 224 ! WRITE(numout,*) ' at_i : ', at_i(12,44) 225 ! WRITE(numout,*) ' v_i : ', v_i(12, 44, 1:jpl) 226 ! WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl) 227 ! !ENDIF 228 229 zhimax = 10._wp 230 ! other categories 231 DO jl = 2, jpl 232 jm = ice_types(jl) 233 DO jj = 1, jpj 234 DO ji = 1, jpi 235 zindb = MAX( rzero, SIGN( rone, ABS(d_a_i_trp(ji,jj,jl)) - epsi10 ) ) 236 ! this correction is very tricky... sometimes, advection gets wrong i don't know why 237 ! it makes problems when the advected volume and concentration do not seem to be 238 ! related with each other 239 ! the new thickness is sometimes very big! 240 ! and sometimes d_a_i_trp and d_v_i_trp have different sign 241 ! which of course is plausible 242 ! but fuck! it fucks everything up :) 243 IF ( (ABS(d_v_i_trp(ji,jj,jl))/MAX(ABS(d_a_i_trp(ji,jj,jl)),epsi10)*zindb .GT. zhimax) & 244 .AND.(v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl),epsi10)*zindb) .GT. zhimax ) THEN 245 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) ) / 2.0 246 a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl) 247 ENDIF 248 END DO ! ji 249 END DO !jj 250 END DO !jl 184 251 185 252 at_i(:,:) = 0._wp … … 187 254 at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 188 255 END DO 256 257 ! !IF( ln_nicep ) THEN 258 ! WRITE(numout,*) ' limupdate1 - after h correction 2 ' 259 ! WRITE(numout,*) ' a_i : ', a_i(12, 44, 1:jpl) 260 ! WRITE(numout,*) ' at_i : ', at_i(12,44) 261 ! WRITE(numout,*) ' v_i : ', v_i(12, 44, 1:jpl) 262 ! WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl) 263 ! !ENDIF 189 264 190 265 !---------------------------------------------------- … … 247 322 !---------------------------------------------------------------------------- 248 323 zindg = tms(ji,jj) * MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) ) ) 249 v_i(ji,jj,jl) = v_i(ji,jj,jl) + zindg * rhosn * v_s(ji,jj,jl) / rau0 250 v_s(ji,jj,jl) = ( 1._wp - zindg ) * v_s(ji,jj,jl) + zindg * v_i(ji,jj,jl) * ( rau0 - rhoic ) / rhosn 324 zdvres = zindg * rhosn * v_s(ji,jj,jl) / rau0 325 v_i(ji,jj,jl) = v_i(ji,jj,jl) + zdvres 326 327 zdvres = zindsn*zindb * ( - zindg * v_s(ji,jj,jl) + zindg * v_i(ji,jj,jl) * ( rau0 - rhoic ) / rhosn ) 328 v_s(ji,jj,jl) = v_s(ji,jj,jl) + zdvres 251 329 252 330 !--- 2.7 Correction to ice concentrations … … 286 364 287 365 !--- 2.13 ice concentration should not exceed amax 366 ! (it should not be the case) 288 367 !----------------------------------------------------- 289 368 DO jj = 1, jpj 290 369 DO ji = 1, jpi 291 370 z_da_ex = MAX( at_i(ji,jj) - amax , 0.0 ) 371 zindb = MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi06 ) ) 292 372 DO jl = 1, jpl 293 zindb = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi04 ) )294 373 z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi06 ) * zindb 295 a_i(ji,jj,jl) = a_i(ji,jj,jl) - z_da_i 296 END DO 297 END DO 298 END DO 374 a_i(ji,jj,jl) = MAX( 0._wp, a_i(ji,jj,jl) - z_da_i ) 375 ! 376 zinda = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) ) 377 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi06 ) * zinda 378 !v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! makes ice shrinken but should not be used 379 END DO 380 END DO 381 END DO 299 382 at_i(:,:) = a_i(:,:,1) 300 383 DO jl = 2, jpl 301 384 at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 302 385 END DO 386 303 387 304 388 ! Final thickness distribution rebinning … … 311 395 ENDIF 312 396 END DO 397 313 398 314 399 !--------------------- … … 356 441 !- check conservation (C Rousset) 357 442 IF (ln_limdiahsb) THEN 443 358 444 zchk_fs = glob_sum( ( fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 359 445 zchk_fw = glob_sum( rdmicif(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b … … 362 448 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic ) 363 449 450 zchk_vmin = glob_min(v_i) 451 zchk_amax = glob_max(SUM(a_i,dim=3)) 452 zchk_amin = glob_min(a_i) 453 364 454 IF(lwp) THEN 365 IF ( 366 IF ( 367 IF ( MINVAL( v_i(:,:,:) ) < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limupdate1) = ',(MINVAL(v_i)* 1.e-3)368 IF ( MAXVAL( SUM(a_i(:,:,:),dim=3) ) > amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax (limupdate1) = ',MAXVAL(SUM(a_i,dim=3))369 IF ( MINVAL( a_i(:,:,:) ) < 0. ) WRITE(numout,*) 'violation a_i<0 (limupdate1) = ',MINVAL(a_i)455 IF ( ABS( zchk_v_i ) > 1.e-5 ) WRITE(numout,*) 'violation volume [m3/day] (limupdate1) = ',(zchk_v_i * 86400.) 456 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limupdate1) = ',(zchk_smv * 86400.) 457 IF ( zchk_vmin < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limupdate1) = ',(zchk_vmin * 1.e-3) 458 IF ( zchk_amax > amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax (limupdate1) = ',zchk_amax 459 IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limupdate1) = ',zchk_amin 370 460 ENDIF 371 461 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
