- 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/limupdate2.F90
r3938 r3963 83 83 INTEGER :: ind_im, layer ! indices for internal melt 84 84 REAL(wp) :: zweight, zesum, zhimax, z_da_i 85 REAL(wp) :: zind b, zindsn, zindic85 REAL(wp) :: zinda, zindb, zindsn, zindic 86 86 REAL(wp) :: zindg, zh, zdvres, zviold2 87 87 REAL(wp) :: zbigvalue, zvsold2, z_da_ex … … 91 91 REAL(wp), POINTER, DIMENSION(:) :: zthick0, zqm0 ! thickness of the layers and heat contents for 92 92 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) 93 93 REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 94 94 ! mass and salt flux (clem) 95 95 REAL(wp), POINTER, DIMENSION(:,:,:) :: zviold, zvsold, zsmvold ! old ice volume... … … 177 177 END DO 178 178 179 180 179 !-------------------------------------- 181 180 ! 2. Review of all pathological cases 182 181 !-------------------------------------- 183 182 !------------------------------------------- 183 ! 2.1) Advection of ice in an ice-free cell 184 !------------------------------------------- 185 ! should be removed since it is treated after dynamics now 186 187 ! !IF( ln_nicep ) THEN 188 ! WRITE(numout,*) ' limupdate2 - before h correction ' 189 ! WRITE(numout,*) ' a_i : ', a_i(12, 44, 1:jpl) 190 ! WRITE(numout,*) ' at_i : ', at_i(12,44) 191 ! WRITE(numout,*) ' v_i : ', v_i(12, 44, 1:jpl) 192 ! WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl) 193 ! !ENDIF 194 195 zhimax = 1._wp 196 ! first category 197 DO jj = 1, jpj 198 DO ji = 1, jpi 199 !--- the thickness of such an ice is often out of bounds 200 !--- thus we recompute a new area while conserving ice volume 201 zat_i_old = SUM( old_a_i(ji,jj,:) ) 202 zindb = MAX( 0._wp, SIGN( 1._wp, ABS( d_a_i_thd(ji,jj,1)) - epsi10 ) ) 203 IF ( ( ABS(d_v_i_thd(ji,jj,1))/MAX(ABS(d_a_i_thd(ji,jj,1)),epsi10)*zindb .GT. zhimax) & 204 .AND.( ( v_i(ji,jj,1)/MAX(a_i(ji,jj,1),epsi10)*zindb) .GT. zhimax ) & 205 .AND.( zat_i_old .LT. 1.e-6 ) ) THEN ! new line 206 ht_i(ji,jj,1) = hi_max(1) / 2.0 207 a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1) 208 ENDIF 209 END DO 210 END DO 211 212 ! !IF( ln_nicep ) THEN 213 ! at_i(:,:) = 0._wp 214 ! DO jl = 1, jpl 215 ! at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 216 ! END DO 217 ! WRITE(numout,*) ' limupdate2 - after h correction 1 ' 218 ! WRITE(numout,*) ' a_i : ', a_i(12, 44, 1:jpl) 219 ! WRITE(numout,*) ' at_i : ', at_i(12,44) 220 ! WRITE(numout,*) ' v_i : ', v_i(12, 44, 1:jpl) 221 ! WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl) 222 ! !ENDIF 223 ! 224 zhimax = 10._wp 225 ! other categories 226 DO jl = 2, jpl 227 jm = ice_types(jl) 228 DO jj = 1, jpj 229 DO ji = 1, jpi 230 zindb = MAX( rzero, SIGN( rone, ABS(d_a_i_thd(ji,jj,jl)) - epsi10 ) ) 231 ! this correction is very tricky... sometimes, advection gets wrong i don't know why 232 ! it makes problems when the advected volume and concentration do not seem to be 233 ! related with each other 234 ! the new thickness is sometimes very big! 235 ! and sometimes d_a_i_trp and d_v_i_trp have different sign 236 ! which of course is plausible 237 ! but fuck! it fucks everything up :) 238 IF ( (ABS(d_v_i_thd(ji,jj,jl))/MAX(ABS(d_a_i_thd(ji,jj,jl)),epsi10)*zindb .GT. zhimax) & 239 .AND.(v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl),epsi10)*zindb) .GT. zhimax ) THEN 240 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 241 a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl) 242 ENDIF 243 END DO ! ji 244 END DO !jj 245 END DO !jl 246 247 at_i(:,:) = 0._wp 248 DO jl = 1, jpl 249 at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 250 END DO 251 252 ! !IF( ln_nicep ) THEN 253 ! WRITE(numout,*) ' limupdate2 - after h correction 2 ' 254 ! WRITE(numout,*) ' a_i : ', a_i(12, 44, 1:jpl) 255 ! WRITE(numout,*) ' at_i : ', at_i(12,44) 256 ! WRITE(numout,*) ' v_i : ', v_i(12, 44, 1:jpl) 257 ! WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl) 258 ! !ENDIF 184 259 185 260 !---------------------------------------------------- … … 363 438 v_i(ji,jj,jl) = v_i(ji,jj,jl) + zdvres 364 439 365 zdvres = - zindg * v_s(ji,jj,jl) + zindg * v_i(ji,jj,jl) * ( rau0 - rhoic ) / rhosn440 zdvres = zindsn*zindb * ( - zindg * v_s(ji,jj,jl) + zindg * v_i(ji,jj,jl) * ( rau0 - rhoic ) / rhosn ) 366 441 v_s(ji,jj,jl) = v_s(ji,jj,jl) + zdvres 367 442 … … 421 496 422 497 !--- 2.13 ice concentration should not exceed amax 498 ! (it should not be the case) 423 499 !----------------------------------------------------- 424 500 DO jj = 1, jpj 425 501 DO ji = 1, jpi 426 502 z_da_ex = MAX( at_i(ji,jj) - amax , 0.0 ) 503 zindb = MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi06 ) ) 427 504 DO jl = 1, jpl 428 zindb = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi04 ) )429 505 z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi06 ) * zindb 430 a_i(ji,jj,jl) = a_i(ji,jj,jl) - z_da_i 506 a_i(ji,jj,jl) = MAX( 0._wp, a_i(ji,jj,jl) - z_da_i ) 507 ! 508 zinda = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) ) 509 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi06 ) * zinda 510 !v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! makes ice shrinken but should not be used 431 511 END DO 432 512 END DO … … 517 597 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic ) 518 598 599 zchk_vmin = glob_min(v_i) 600 zchk_amax = glob_max(SUM(a_i,dim=3)) 601 zchk_amin = glob_min(a_i) 602 519 603 IF(lwp) THEN 520 IF ( 521 IF ( 522 IF ( MINVAL( v_i(:,:,:) ) < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limupdate2) = ',(MINVAL(v_i)* 1.e-3)523 IF ( MAXVAL( SUM(a_i(:,:,:),dim=3) ) > amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax (limupdate2) = ',MAXVAL(SUM(a_i,dim=3))524 IF ( MINVAL( a_i(:,:,:) ) < 0. ) WRITE(numout,*) 'violation a_i<0 (limupdate2) = ',MINVAL(a_i)604 IF ( ABS( zchk_v_i ) > 1.e-5 ) WRITE(numout,*) 'violation volume [m3/day] (limupdate2) = ',(zchk_v_i * 86400.) 605 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limupdate2) = ',(zchk_smv * 86400.) 606 IF ( zchk_vmin < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limupdate2) = ',(zchk_vmin * 1.e-3) 607 IF ( zchk_amax > amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax (limupdate2) = ',zchk_amax 608 IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limupdate2) = ',zchk_amin 525 609 ENDIF 526 610 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
