- Timestamp:
- 2013-12-11T15:38:42+01:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90
r4099 r4332 48 48 PUBLIC lim_update1 ! routine called by ice_step 49 49 50 REAL(wp) :: epsi06 = 1.e-06_wp ! module constants51 REAL(wp) :: epsi04 = 1.e-04_wp ! - -52 REAL(wp) :: epsi03 = 1.e-03_wp ! - -53 50 REAL(wp) :: epsi10 = 1.e-10_wp ! - - 54 REAL(wp) :: epsi16 = 1.e-16_wp ! - -55 REAL(wp) :: epsi20 = 1.e-20_wp ! - -56 51 REAL(wp) :: rzero = 0._wp ! - - 57 52 REAL(wp) :: rone = 1._wp ! - - … … 108 103 !------------------------------------------------------------------------------ 109 104 110 !--------------------- 111 ! Ice dynamics 112 !--------------------- 113 u_ice(:,:) = u_ice(:,:) + d_u_ice_dyn(:,:) 114 v_ice(:,:) = v_ice(:,:) + d_v_ice_dyn(:,:) 115 116 !----------------------------- 117 ! Update ice and snow volumes 118 !----------------------------- 119 DO jl = 1, jpl 120 v_i(:,:,jl) = v_i(:,:,jl) + d_v_i_trp(:,:,jl) 121 v_s(:,:,jl) = v_s(:,:,jl) + d_v_s_trp(:,:,jl) 122 END DO 123 124 125 !--------------------------------------------- 126 ! Ice concentration and ice heat content 127 !--------------------------------------------- 128 a_i (:,:,:) = a_i (:,:,:) + d_a_i_trp(:,:,:) 129 e_i(:,:,:,:) = e_i(:,:,:,:) + d_e_i_trp(:,:,:,:) 130 131 !------------------------------ 132 ! Snow temperature and ice age 133 !------------------------------ 134 e_s (:,:,:,:) = e_s (:,:,:,:) + d_e_s_trp (:,:,:,:) 135 oa_i(:,:,:) = oa_i(:,:,:) + d_oa_i_trp(:,:,:) 136 137 !-------------- 138 ! Ice salinity 139 !-------------- 140 IF( num_sal == 2 ) THEN ! general case 141 smv_i(:,:,:) = smv_i(:,:,:) + d_smv_i_trp(:,:,:) 142 ENDIF 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(:,:,:) 143 118 144 119 ! mass and salt flux init (clem) … … 160 135 CALL lim_var_glo2eqv 161 136 162 !---------------------------------163 ! Classify the pathological cases164 !---------------------------------165 ! (1) v_i (new) > 0; d_v_i_thd + v_i(old) > 0 (easy case)166 ! (3) v_i (new) < 0; d_v_i_thd + v_i(old) > 0 (combined total ablation)167 ! (5) v_i (old) = 0; d_v_i_trp > 0 (advection of ice in a free-cell)168 !169 DO jl = 1, jpl170 DO jj = 1, jpj171 DO ji = 1, jpi172 patho_case(ji,jj,jl) = 1173 IF( v_i(ji,jj,jl) .LT. 0.0 ) THEN174 patho_case(ji,jj,jl) = 3175 ENDIF176 IF( ( old_v_i(ji,jj,jl) .LE. epsi10 ) .AND. ( d_v_i_trp(ji,jj,jl) .GT. epsi06 ) ) THEN177 patho_case(ji,jj,jl) = 5 ! advection of ice in an ice-free cell178 ENDIF179 END DO180 END DO181 END DO182 183 137 !-------------------------------------- 184 138 ! 2. Review of all pathological cases 185 139 !-------------------------------------- 186 140 141 ! clem: useless now 187 142 !------------------------------------------- 188 143 ! 2.1) Advection of ice in an ice-free cell 189 144 !------------------------------------------- 190 145 ! should be removed since it is treated after dynamics now 191 192 ! !IF( ln_nicep ) THEN 193 ! WRITE(numout,*) ' limupdate1 - before h correction ' 194 ! WRITE(numout,*) ' a_i : ', a_i(12, 44, 1:jpl) 195 ! WRITE(numout,*) ' at_i : ', at_i(12,44) 196 ! WRITE(numout,*) ' v_i : ', v_i(12, 44, 1:jpl) 197 ! WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl) 198 ! !ENDIF 146 ! zhimax = 5._wp 147 ! ! first category 148 ! DO jj = 1, jpj 149 ! DO ji = 1, jpi 150 ! !--- the thickness of such an ice is often out of bounds 151 ! !--- thus we recompute a new area while conserving ice volume 152 ! 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 line 157 ! ht_i(ji,jj,1) = hi_max(1) * 0.5_wp 158 ! a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1) 159 ! ENDIF 160 ! END DO 161 ! END DO 199 162 ! 200 zhimax = .3_wp 201 ! first category 202 DO jj = 1, jpj 203 DO ji = 1, jpi 204 !--- the thickness of such an ice is often out of bounds 205 !--- thus we recompute a new area while conserving ice volume 206 zat_i_old = SUM( old_a_i(ji,jj,:) ) 207 zindb = MAX( 0._wp, SIGN( 1._wp, ABS( d_a_i_trp(ji,jj,1)) - epsi10 ) ) 208 IF ( ( ABS(d_v_i_trp(ji,jj,1))/MAX(ABS(d_a_i_trp(ji,jj,1)),epsi10)*zindb .GT. zhimax) & 209 .AND.( ( v_i(ji,jj,1)/MAX(a_i(ji,jj,1),epsi10)*zindb) .GT. zhimax ) & 210 .AND.( zat_i_old .LT. 1.e-6 ) ) THEN ! new line 211 ht_i(ji,jj,1) = hi_max(1) / 2.0 212 a_i (ji,jj,1) = v_i(ji,jj,1) / ht_i(ji,jj,1) 213 ENDIF 214 END DO 215 END DO 216 217 218 ! !IF( ln_nicep ) THEN 219 ! at_i(:,:) = 0._wp 220 ! DO jl = 1, jpl 221 ! at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 222 ! END DO 223 ! WRITE(numout,*) ' limupdate1 - after h correction 1 ' 224 ! WRITE(numout,*) ' a_i : ', a_i(12, 44, 1:jpl) 225 ! WRITE(numout,*) ' at_i : ', at_i(12,44) 226 ! WRITE(numout,*) ' v_i : ', v_i(12, 44, 1:jpl) 227 ! WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl) 228 ! !ENDIF 229 230 zhimax = 1._wp 231 ! other categories 232 DO jl = 2, jpl 233 jm = ice_types(jl) 234 DO jj = 1, jpj 235 DO ji = 1, jpi 236 zindb = MAX( rzero, SIGN( rone, ABS(d_a_i_trp(ji,jj,jl)) - epsi10 ) ) 237 ! this correction is very tricky... sometimes, advection gets wrong i don't know why 238 ! it makes problems when the advected volume and concentration do not seem to be 239 ! related with each other 240 ! the new thickness is sometimes very big! 241 ! and sometimes d_a_i_trp and d_v_i_trp have different sign 242 ! which of course is plausible 243 ! but fuck! it fucks everything up :) 244 IF ( (ABS(d_v_i_trp(ji,jj,jl))/MAX(ABS(d_a_i_trp(ji,jj,jl)),epsi10)*zindb .GT. zhimax) & 245 .AND.(v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl),epsi10)*zindb) .GT. zhimax ) THEN 246 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 247 a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl) 248 ENDIF 249 END DO ! ji 250 END DO !jj 251 END DO !jl 163 ! zhimax = 20._wp 164 ! ! other categories 165 ! DO jl = 2, jpl 166 ! jm = ice_types(jl) 167 ! DO jj = 1, jpj 168 ! DO ji = 1, jpi 169 ! 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 why 171 ! ! it makes problems when the advected volume and concentration do not seem to be 172 ! ! related with each other 173 ! ! the new thickness is sometimes very big! 174 ! ! and sometimes d_a_i_trp and d_v_i_trp have different sign 175 ! ! which of course is plausible 176 ! ! 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 ) THEN 179 ! 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 180 ! a_i (ji,jj,jl) = v_i(ji,jj,jl) / ht_i(ji,jj,jl) 181 ! ENDIF 182 ! END DO ! ji 183 ! END DO !jj 184 ! END DO !jl 252 185 253 186 at_i(:,:) = 0._wp … … 255 188 at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 256 189 END DO 257 258 ! !IF( ln_nicep ) THEN259 ! WRITE(numout,*) ' limupdate1 - after h correction 2 '260 ! WRITE(numout,*) ' a_i : ', a_i(12, 44, 1:jpl)261 ! WRITE(numout,*) ' at_i : ', at_i(12,44)262 ! WRITE(numout,*) ' v_i : ', v_i(12, 44, 1:jpl)263 ! WRITE(numout,*) ' ht_i : ', ht_i(12, 44, 1:jpl)264 ! !ENDIF265 190 266 191 !---------------------------------------------------- … … 285 210 286 211 !switches 287 zindb = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi 06) )212 zindb = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi10 ) ) 288 213 !switch = 1 if a_i > 1e-06 and 0 if not 289 zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi 06 ) ) !=1 if hs > 1e-6and 0 if not290 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi 04 ) ) !=1 if hi > 1e-3and 0 if not214 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 291 216 ! bug fix 25 avril 2007 292 217 zindb = zindb*zindic … … 332 257 !-------------------------------------------- 333 258 ! if greater than 0, ice concentration cannot be smaller than 1e-10 334 !clem a_i(ji,jj,jl) = zindb * MAX(zindsn, zindic) * MAX( a_i(ji,jj,jl), epsi06 )335 259 a_i(ji,jj,jl) = zindb * a_i(ji,jj,jl) 336 260 … … 352 276 DO jj = 1, jpj 353 277 DO ji = 1, jpi 354 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi 04) )278 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) 355 279 e_i(ji,jj,jk,jl)= zindic * ( MIN ( MAX ( 0.0, e_i(ji,jj,jk,jl) ), zbigvalue ) ) 356 280 END DO ! ji … … 370 294 DO ji = 1, jpi 371 295 z_da_ex = MAX( at_i(ji,jj) - amax , 0.0 ) 372 zindb = MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi 06) )296 zindb = MAX( rzero, SIGN( rone, at_i(ji,jj) - epsi10 ) ) 373 297 DO jl = 1, jpl 374 z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi 06) * zindb298 z_da_i = a_i(ji,jj,jl) * z_da_ex / MAX( at_i(ji,jj), epsi10 ) * zindb 375 299 a_i(ji,jj,jl) = MAX( 0._wp, a_i(ji,jj,jl) - z_da_i ) 376 300 ! 377 zinda = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi 06) )378 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi 06) * zinda301 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 379 303 !v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl) ! makes ice shrinken but should not be used 380 304 END DO … … 412 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) ) 413 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) ) 414 i_ice_switch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) + epsi20) )415 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)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) 416 340 END DO ! ji 417 341 END DO ! jj … … 512 436 CALL prt_ctl(tab2d_1=old_oa_i (:,:,jl) , clinfo1= ' lim_update1 : old_oa_i : ') 513 437 CALL prt_ctl(tab2d_1=d_oa_i_trp (:,:,jl) , clinfo1= ' lim_update1 : d_oa_i_trp : ') 514 CALL prt_ctl(tab2d_1=REAL(patho_case(:,:,jl)) , clinfo1= ' lim_update1 : Path. case : ')515 438 516 439 DO jk = 1, nlay_i
Note: See TracChangeset
for help on using the changeset viewer.
