- Timestamp:
- 2013-12-11T18:34:00+01:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
r4293 r4333 45 45 PUBLIC lim_update2 ! routine called by ice_step 46 46 47 REAL(wp) :: epsi06 = 1.e-06_wp ! module constants48 REAL(wp) :: epsi04 = 1.e-04_wp ! - -49 47 REAL(wp) :: epsi10 = 1.e-10_wp ! - - 50 REAL(wp) :: epsi16 = 1.e-16_wp ! - -51 REAL(wp) :: epsi20 = 1.e-20_wp ! - -52 48 REAL(wp) :: rzero = 0._wp ! - - 53 49 REAL(wp) :: rone = 1._wp ! - - … … 102 98 CALL wrk_alloc( jpi,jpj,jpl,zviold, zvsold, zsmvold ) ! clem 103 99 104 !------------------------------------------------------------------------------ 105 ! 1. Update of Global variables | 106 !------------------------------------------------------------------------------ 107 108 !----------------------------- 109 ! Update ice and snow volumes 110 !----------------------------- 111 DO jl = 1, jpl 112 v_i(:,:,jl) = v_i(:,:,jl) + d_v_i_thd(:,:,jl) 113 v_s(:,:,jl) = v_s(:,:,jl) + d_v_s_thd(:,:,jl) 114 END DO 115 116 !--------------------------------------------- 117 ! Ice concentration and ice heat content 118 !--------------------------------------------- 119 a_i (:,:,:) = a_i (:,:,:) + d_a_i_thd(:,:,:) 120 e_i(:,:,:,:) = e_i(:,:,:,:) + d_e_i_thd(:,:,:,:) 121 122 !------------------------------ 123 ! Snow temperature and ice age 124 !------------------------------ 125 e_s (:,:,:,:) = e_s (:,:,:,:) + d_e_s_thd (:,:,:,:) 126 oa_i(:,:,:) = oa_i(:,:,:) + d_oa_i_thd(:,:,:) 127 128 !-------------- 129 ! Ice salinity 130 !-------------- 131 IF( num_sal == 2 ) THEN ! general case 132 smv_i(:,:,:) = smv_i(:,:,:) + d_smv_i_thd(:,:,:) 133 ENDIF 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 134 114 135 115 ! mass and salt flux init (clem) … … 151 131 CALL lim_var_glo2eqv 152 132 153 !---------------------------------154 ! Classify the pathological cases155 !---------------------------------156 ! (1) v_i (new) > 0; d_v_i_thd + v_i(old) > 0 (easy case)157 ! (2) v_i (new) > 0; d_v_i_thd + v_i(old) = 0 (total thermodynamic ablation)158 ! (3) v_i (new) < 0; d_v_i_thd + v_i(old) > 0 (combined total ablation)159 ! (4) v_i (new) < 0; d_v_i_thd + v_i(old) = 0 (total thermodynamic ablation160 ! with negative advection, very pathological )161 !162 DO jl = 1, jpl163 DO jj = 1, jpj164 DO ji = 1, jpi165 patho_case(ji,jj,jl) = 1166 IF( v_i(ji,jj,jl) .GE. 0.0 ) THEN167 IF ( old_v_i(ji,jj,jl) + d_v_i_thd(ji,jj,jl) .LT. epsi10 ) THEN168 patho_case(ji,jj,jl) = 2169 ENDIF170 ELSE171 patho_case(ji,jj,jl) = 3172 IF( old_v_i(ji,jj,jl) + d_v_i_thd(ji,jj,jl) .LT. epsi10 ) THEN173 patho_case(ji,jj,jl) = 4174 ENDIF175 ENDIF176 END DO177 END DO178 END DO179 180 133 !-------------------------------------- 181 134 ! 2. Review of all pathological cases 182 135 !-------------------------------------- 136 137 ! clem: useless now 183 138 !------------------------------------------- 184 139 ! 2.1) Advection of ice in an ice-free cell 185 140 !------------------------------------------- 186 141 ! should be removed since it is treated after dynamics now 187 188 ! !IF( ln_nicep ) THEN 189 ! WRITE(numout,*) ' limupdate2 - before h correction ' 190 ! WRITE(numout,*) ' a_i : ', a_i(12, 44, 1:jpl) 191 ! WRITE(numout,*) ' at_i : ', at_i(12,44) 192 ! WRITE(numout,*) ' v_i : ', v_i(12, 44, 1:jpl) 193 ! WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl) 194 ! !ENDIF 195 196 zhimax = 5._wp 197 ! first category 198 DO jj = 1, jpj 199 DO ji = 1, jpi 200 !--- the thickness of such an ice is often out of bounds 201 !--- thus we recompute a new area while conserving ice volume 202 zat_i_old = SUM( old_a_i(ji,jj,:) ) 203 zindb = MAX( 0._wp, SIGN( 1._wp, ABS( d_a_i_thd(ji,jj,1) ) - epsi10 ) ) 204 IF ( ( ABS( d_v_i_thd(ji,jj,1) ) / MAX( ABS( d_a_i_thd(ji,jj,1) ),epsi10 ) * zindb .GT. zhimax ) & 205 & .AND. ( ( v_i(ji,jj,1) / MAX( a_i(ji,jj,1), epsi10 ) * zindb ) .GT. zhimax ) & 206 & .AND. ( zat_i_old .LT. 1.e-6 ) ) THEN ! new line 207 ht_i(ji,jj,1) = hi_max(1) * 0.5_wp 208 a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1) 209 ENDIF 210 END DO 211 END DO 212 213 ! !IF( ln_nicep ) THEN 214 ! at_i(:,:) = 0._wp 215 ! DO jl = 1, jpl 216 ! at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 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 217 156 ! END DO 218 ! WRITE(numout,*) ' limupdate2 - after h correction 1 ' 219 ! WRITE(numout,*) ' a_i : ', a_i(12, 44, 1:jpl) 220 ! WRITE(numout,*) ' at_i : ', at_i(12,44) 221 ! WRITE(numout,*) ' v_i : ', v_i(12, 44, 1:jpl) 222 ! WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl) 223 ! !ENDIF 224 ! 225 zhimax = 20._wp 226 ! other categories 227 DO jl = 2, jpl 228 jm = ice_types(jl) 229 DO jj = 1, jpj 230 DO ji = 1, jpi 231 zindb = MAX( rzero, SIGN( rone, ABS( d_a_i_thd(ji,jj,jl)) - epsi10 ) ) 232 ! this correction is very tricky... sometimes, advection gets wrong i don't know why 233 ! it makes problems when the advected volume and concentration do not seem to be 234 ! related with each other 235 ! the new thickness is sometimes very big! 236 ! and sometimes d_a_i_trp and d_v_i_trp have different sign 237 ! which of course is plausible 238 ! but fuck! it fucks everything up :) 239 IF ( ( ABS( d_v_i_thd(ji,jj,jl) ) / MAX( ABS( d_a_i_thd(ji,jj,jl) ), epsi10 ) * zindb .GT. zhimax ) & 240 & .AND. ( v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zindb ) .GT. zhimax ) THEN 241 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 242 a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl) 243 ENDIF 244 END DO ! ji 245 END DO !jj 246 END DO !jl 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 247 181 248 182 at_i(:,:) = 0._wp … … 250 184 at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 251 185 END DO 252 253 ! !IF( ln_nicep ) THEN254 ! WRITE(numout,*) ' limupdate2 - after h correction 2 '255 ! WRITE(numout,*) ' a_i : ', a_i(12, 44, 1:jpl)256 ! WRITE(numout,*) ' at_i : ', at_i(12,44)257 ! WRITE(numout,*) ' v_i : ', v_i(12, 44, 1:jpl)258 ! WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl)259 ! !ENDIF260 186 261 187 !---------------------------------------------------- … … 304 230 zesum = zesum + e_i(ji,jj,jk,jl) 305 231 END DO 306 !clem IF (ind_im .LT. nlay_i ) THEN307 !clem smv_i(ji,jj,jl) = smv_i(ji,jj,jl) / ht_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - REAL(ind_im)*ht_i(ji,jj,jl) / REAL(nlay_i) )308 !clem ENDIF309 232 ht_i(ji,jj,jl) = ht_i(ji,jj,jl) - REAL(ind_im)*ht_i(ji,jj,jl) / REAL(nlay_i) 310 233 v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) … … 368 291 DO ji = 1, jpi 369 292 ! snow energy of melting 370 ze_s = e_s(ji,jj,1,jl) * unit_fac / area(ji,jj) / MAX( v_s(ji,jj,jl), 1.0e-6 ) ! 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 371 295 372 296 ! If snow energy of melting smaller then Lf … … 401 325 402 326 !switches 403 zindb = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi 06) )327 zindb = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi10 ) ) 404 328 !switch = 1 if a_i > 1e-06 and 0 if not 405 zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi 06 ) ) !=1 if hs > 1e-6and 0 if not406 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi 04 ) ) !=1 if hi > 1e-3and 0 if not329 zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) ) !=1 if hs > 1e-10 and 0 if not 330 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) !=1 if hi > 1e-10 and 0 if not 407 331 ! bug fix 25 avril 2007 408 332 zindb = zindb*zindic … … 444 368 !--- 2.7 Correction to ice concentrations 445 369 !-------------------------------------------- 446 !clem a_i(ji,jj,jl) = zindb * MAX(zindsn, zindic) * MAX( a_i(ji,jj,jl), epsi06 )447 370 a_i(ji,jj,jl) = zindb * a_i(ji,jj,jl) 448 371 … … 464 387 DO jj = 1, jpj 465 388 DO ji = 1, jpi 466 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi 04) )389 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) 467 390 e_i(ji,jj,jk,jl)= zindic * ( MIN ( MAX ( 0.0, e_i(ji,jj,jk,jl) ), zbigvalue ) ) 468 391 END DO ! ji … … 478 401 !--- 2.12 Constrain the thickness of the smallest category above 5 cm 479 402 !---------------------------------------------------------------------- 480 zindb = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi 06) )481 ht_i(ji,jj,jl) = zindb*v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl), epsi 06)482 zh = MAX( rone , zindb * hiclim / MAX( ht_i(ji,jj,jl) , epsi 20 ) )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 ) ) 483 406 ht_s(ji,jj,jl) = ht_s(ji,jj,jl)* zh 484 407 ht_i(ji,jj,jl) = ht_i(ji,jj,jl)* zh … … 502 425 DO ji = 1, jpi 503 426 z_da_ex = MAX( at_i(ji,jj) - amax , 0.0 ) 504 zindb = MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi 06) )427 zindb = MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi10 ) ) 505 428 DO jl = 1, jpl 506 z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi 06) * zindb429 z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi10 ) * zindb 507 430 a_i(ji,jj,jl) = MAX( 0._wp, a_i(ji,jj,jl) - z_da_i ) 508 431 ! 509 zinda = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi 06) )510 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi 06) * zinda432 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 511 434 !v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! makes ice shrinken but should not be used 512 435 END DO … … 542 465 !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) ) 543 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) ) 544 i_ice_switch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) + epsi20) )545 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)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) 546 469 END DO ! ji 547 470 END DO ! jj … … 663 586 CALL prt_ctl(tab2d_1=old_oa_i (:,:,jl) , clinfo1= ' lim_update2 : old_oa_i : ') 664 587 CALL prt_ctl(tab2d_1=d_oa_i_thd (:,:,jl) , clinfo1= ' lim_update2 : d_oa_i_thd : ') 665 CALL prt_ctl(tab2d_1=REAL(patho_case(:,:,jl)) , clinfo1= ' lim_update2 : Path. case : ')666 588 667 589 DO jk = 1, nlay_i
Note: See TracChangeset
for help on using the changeset viewer.
