Changeset 5682 for branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
- Timestamp:
- 2015-08-12T17:46:45+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
r4990 r5682 13 13 !! 'key_lim3' : LIM3 sea-ice model 14 14 !!---------------------------------------------------------------------- 15 !! lim_itd_th : thermodynamics of ice thickness distribution16 15 !! lim_itd_th_rem : 17 16 !! lim_itd_th_reb : … … 25 24 USE thd_ice ! LIM-3 thermodynamic variables 26 25 USE ice ! LIM-3 variables 27 USE par_ice ! LIM-3 parameters28 USE limthd_lac ! LIM-3 lateral accretion29 26 USE limvar ! LIM-3 variables 30 USE limcons ! LIM-3 conservation31 27 USE prtctl ! Print control 32 28 USE in_out_manager ! I/O manager … … 34 30 USE wrk_nemo ! work arrays 35 31 USE lib_fortran ! to use key_nosignedzero 36 USE timing ! Timing 37 USE limcons ! conservation tests 32 USE limcons ! conservation tests 38 33 39 34 IMPLICIT NONE 40 35 PRIVATE 41 36 42 PUBLIC lim_itd_th ! called by ice_stp43 37 PUBLIC lim_itd_th_rem 44 38 PUBLIC lim_itd_th_reb 45 PUBLIC lim_itd_fitline46 PUBLIC lim_itd_shiftice47 39 48 40 !!---------------------------------------------------------------------- … … 53 45 CONTAINS 54 46 55 SUBROUTINE lim_itd_th( kt )56 !!------------------------------------------------------------------57 !! *** ROUTINE lim_itd_th ***58 !!59 !! ** Purpose : computes the thermodynamics of ice thickness distribution60 !!61 !! ** Method :62 !!------------------------------------------------------------------63 INTEGER, INTENT(in) :: kt ! time step index64 !65 INTEGER :: ji, jj, jk, jl ! dummy loop index66 !67 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b68 !!------------------------------------------------------------------69 IF( nn_timing == 1 ) CALL timing_start('limitd_th')70 71 ! conservation test72 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_th', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)73 74 IF( kt == nit000 .AND. lwp ) THEN75 WRITE(numout,*)76 WRITE(numout,*) 'lim_itd_th : Thermodynamics of the ice thickness distribution'77 WRITE(numout,*) '~~~~~~~~~~~'78 ENDIF79 80 !------------------------------------------------------------------------------|81 ! 1) Transport of ice between thickness categories. |82 !------------------------------------------------------------------------------|83 ! Given thermodynamic growth rates, transport ice between84 ! thickness categories.85 IF( jpl > 1 ) CALL lim_itd_th_rem( 1, jpl, kt )86 !87 CALL lim_var_glo2eqv ! only for info88 CALL lim_var_agg(1)89 90 !------------------------------------------------------------------------------|91 ! 3) Add frazil ice growing in leads.92 !------------------------------------------------------------------------------|93 CALL lim_thd_lac94 CALL lim_var_glo2eqv ! only for info95 96 IF(ln_ctl) THEN ! Control print97 CALL prt_ctl_info(' ')98 CALL prt_ctl_info(' - Cell values : ')99 CALL prt_ctl_info(' ~~~~~~~~~~~~~ ')100 CALL prt_ctl(tab2d_1=area , clinfo1=' lim_itd_th : cell area :')101 CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_th : at_i :')102 CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_th : vt_i :')103 CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_th : vt_s :')104 DO jl = 1, jpl105 CALL prt_ctl_info(' ')106 CALL prt_ctl_info(' - Category : ', ivar1=jl)107 CALL prt_ctl_info(' ~~~~~~~~~~')108 CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_itd_th : a_i : ')109 CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_itd_th : ht_i : ')110 CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_itd_th : ht_s : ')111 CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_itd_th : v_i : ')112 CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_itd_th : v_s : ')113 CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_itd_th : e_s : ')114 CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_itd_th : t_su : ')115 CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_itd_th : t_snow : ')116 CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_itd_th : sm_i : ')117 CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_itd_th : smv_i : ')118 DO jk = 1, nlay_i119 CALL prt_ctl_info(' ')120 CALL prt_ctl_info(' - Layer : ', ivar1=jk)121 CALL prt_ctl_info(' ~~~~~~~')122 CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_th : t_i : ')123 CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_th : e_i : ')124 END DO125 END DO126 ENDIF127 !128 ! conservation test129 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_th', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)130 !131 IF( nn_timing == 1 ) CALL timing_stop('limitd_th')132 END SUBROUTINE lim_itd_th133 !134 135 47 SUBROUTINE lim_itd_th_rem( klbnd, kubnd, kt ) 136 48 !!------------------------------------------------------------------ … … 153 65 REAL(wp) :: zx1, zwk1, zdh0, zetamin, zdamax ! local scalars 154 66 REAL(wp) :: zx2, zwk2, zda0, zetamax ! - - 155 REAL(wp) :: zx3 , zareamin ! - -67 REAL(wp) :: zx3 156 68 CHARACTER (len = 15) :: fieldid 157 69 … … 179 91 !!------------------------------------------------------------------ 180 92 181 CALL wrk_alloc( jpi,jpj, zremap_flag ) ! integer182 CALL wrk_alloc( jpi,jpj,jpl-1, zdonor ) ! integer93 CALL wrk_alloc( jpi,jpj, zremap_flag ) 94 CALL wrk_alloc( jpi,jpj,jpl-1, zdonor ) 183 95 CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es ) 184 96 CALL wrk_alloc( jpi,jpj,jpl-1, zdaice, zdvice ) 185 97 CALL wrk_alloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 ) 186 98 CALL wrk_alloc( (jpi+1)*(jpj+1), zvetamin, zvetamax ) 187 CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) ! integer99 CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 188 100 CALL wrk_alloc( jpi,jpj, zhb0,zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final ) 189 190 zareamin = epsi10 !minimum area in thickness categories tolerated by the conceptors of the model191 101 192 102 !!---------------------------------------------------------------------------------------------- … … 216 126 DO jj = 1, jpj 217 127 DO ji = 1, jpi 218 rswitch = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i(ji,jj,jl) +epsi10 ) ) !0 if no ice and 1 if yes128 rswitch = MAX( 0.0, SIGN( 1.0, a_i(ji,jj,jl) - epsi10 ) ) !0 if no ice and 1 if yes 219 129 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * rswitch 220 rswitch = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i_b(ji,jj,jl) + epsi10) ) !0 if no ice and 1 if yes130 rswitch = MAX( 0.0, SIGN( 1.0, a_i_b(ji,jj,jl) - epsi10) ) 221 131 zht_i_b(ji,jj,jl) = v_i_b(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) * rswitch 222 IF( a_i(ji,jj,jl) > epsi10 ) zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl) 132 IF( a_i(ji,jj,jl) > epsi10 ) zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl) ! clem: useless IF statement? 223 133 END DO 224 134 END DO … … 239 149 DO jj = 1, jpj 240 150 DO ji = 1, jpi 241 IF ( at_i(ji,jj) .gt. zareamin) THEN151 IF ( at_i(ji,jj) > epsi10 ) THEN 242 152 nbrem = nbrem + 1 243 153 nind_i(nbrem) = ji … … 247 157 zremap_flag(ji,jj) = 0 248 158 ENDIF 249 END DO !ji250 END DO !jj159 END DO 160 END DO 251 161 252 162 !----------------------------------------------------------------------------------------------- … … 254 164 !----------------------------------------------------------------------------------------------- 255 165 !- 4.1 Compute category boundaries 256 ! Tricky trick see limitd_me.F90257 ! will be soon removed, CT258 ! hi_max(kubnd) = 99.259 166 zhbnew(:,:,:) = 0._wp 260 167 … … 265 172 ! 266 173 zhbnew(ii,ij,jl) = hi_max(jl) 267 IF ( a_i_b(ii,ij,jl) > epsi10 .AND. a_i_b(ii,ij,jl+1) > epsi10 ) THEN174 IF ( a_i_b(ii,ij,jl) > epsi10 .AND. a_i_b(ii,ij,jl+1) > epsi10 ) THEN 268 175 !interpolate between adjacent category growth rates 269 176 zslope = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_b(ii,ij,jl+1) - zht_i_b(ii,ij,jl) ) 270 177 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_b(ii,ij,jl) ) 271 ELSEIF 178 ELSEIF( a_i_b(ii,ij,jl) > epsi10) THEN 272 179 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) 273 ELSEIF 180 ELSEIF( a_i_b(ii,ij,jl+1) > epsi10) THEN 274 181 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1) 275 182 ENDIF … … 280 187 ii = nind_i(ji) 281 188 ij = nind_j(ji) 282 IF( a_i(ii,ij,jl) > epsi10 .AND. ht_i(ii,ij,jl) >= zhbnew(ii,ij,jl) ) THEN 189 190 ! clem: we do not want ht_i to be too close to either HR or HL otherwise a division by nearly 0 is possible 191 ! in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 192 IF ( a_i(ii,ij,jl ) > epsi10 .AND. ht_i(ii,ij,jl ) > ( zhbnew(ii,ij,jl) - epsi10 ) ) THEN 283 193 zremap_flag(ii,ij) = 0 284 ELSEIF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) < = zhbnew(ii,ij,jl) ) THEN194 ELSEIF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) < ( zhbnew(ii,ij,jl) + epsi10 ) ) THEN 285 195 zremap_flag(ii,ij) = 0 286 196 ENDIF 287 197 288 198 !- 4.3 Check that each zhbnew does not exceed maximal values hi_max 199 IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0 289 200 IF( zhbnew(ii,ij,jl) > hi_max(jl+1) ) zremap_flag(ii,ij) = 0 290 IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0 291 END DO 292 293 END DO !jl 201 ! clem bug: why is not the following instead? 202 !!IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0 203 !!IF( zhbnew(ii,ij,jl) > hi_max(jl ) ) zremap_flag(ii,ij) = 0 204 205 END DO 206 207 END DO 294 208 295 209 !----------------------------------------------------------------------------------------------- … … 312 226 DO jj = 1, jpj 313 227 DO ji = 1, jpi 314 zhb0(ji,jj) = hi_max(0) ! 0eme 315 zhb1(ji,jj) = hi_max(1) ! 1er 316 317 zhbnew(ji,jj,klbnd-1) = 0._wp 228 zhb0(ji,jj) = hi_max(0) 229 zhb1(ji,jj) = hi_max(1) 318 230 319 231 IF( a_i(ji,jj,kubnd) > epsi10 ) THEN 320 zhbnew(ji,jj,kubnd) = 3._wp * ht_i(ji,jj,kubnd) - 2._wp * zhbnew(ji,jj,kubnd-1)232 zhbnew(ji,jj,kubnd) = MAX( hi_max(kubnd-1), 3._wp * ht_i(ji,jj,kubnd) - 2._wp * zhbnew(ji,jj,kubnd-1) ) 321 233 ELSE 322 zhbnew(ji,jj,kubnd) = hi_max(kubnd) 323 !!? clem bug: since hi_max(jpl)=99, this limit is very high 324 !!? but I think it is erased in fitline subroutine 325 ENDIF 326 327 IF( zhbnew(ji,jj,kubnd) < hi_max(kubnd-1) ) zhbnew(ji,jj,kubnd) = hi_max(kubnd-1) 328 329 END DO !jj 330 END DO !jj 234 !clem bug zhbnew(ji,jj,kubnd) = hi_max(kubnd) 235 zhbnew(ji,jj,kubnd) = hi_max(kubnd-1) ! not used anyway 236 ENDIF 237 238 ! clem: we do not want ht_i_b to be too close to either HR or HL otherwise a division by nearly 0 is possible 239 ! in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 240 IF ( zht_i_b(ji,jj,klbnd) < ( zhb0(ji,jj) + epsi10 ) ) THEN 241 zremap_flag(ji,jj) = 0 242 ELSEIF( zht_i_b(ji,jj,klbnd) > ( zhb1(ji,jj) - epsi10 ) ) THEN 243 zremap_flag(ji,jj) = 0 244 ENDIF 245 246 END DO 247 END DO 331 248 332 249 !----------------------------------------------------------------------------------------------- … … 334 251 !----------------------------------------------------------------------------------------------- 335 252 !- 7.1 g(h) for category 1 at start of time step 336 CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd), & 337 & g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd), & 253 CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd), g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd), & 338 254 & hR(:,:,klbnd), zremap_flag ) 339 255 … … 343 259 ij = nind_j(ji) 344 260 345 !ji346 IF (a_i(ii,ij,klbnd) .gt. epsi10) THEN 261 IF( a_i(ii,ij,klbnd) > epsi10 ) THEN 262 347 263 zdh0 = zdhice(ii,ij,klbnd) !decrease of ice thickness in the lower category 348 ! ji, a_i > epsi10 349 IF (zdh0 .lt. 0.0) THEN !remove area from category 1 350 ! ji, a_i > epsi10; zdh0 < 0 351 zdh0 = MIN(-zdh0,hi_max(klbnd)) 352 264 265 IF( zdh0 < 0.0 ) THEN !remove area from category 1 266 zdh0 = MIN( -zdh0, hi_max(klbnd) ) 353 267 !Integrate g(1) from 0 to dh0 to estimate area melted 354 zetamax = MIN(zdh0,hR(ii,ij,klbnd)) - hL(ii,ij,klbnd) 355 IF (zetamax.gt.0.0) THEN 356 zx1 = zetamax 357 zx2 = 0.5 * zetamax*zetamax 358 zda0 = g1(ii,ij,klbnd) * zx2 + g0(ii,ij,klbnd) * zx1 !ice area removed 359 ! Constrain new thickness <= ht_i 360 zdamax = a_i(ii,ij,klbnd) * & 361 (1.0 - ht_i(ii,ij,klbnd)/zht_i_b(ii,ij,klbnd)) ! zdamax > 0 362 !ice area lost due to melting of thin ice 363 zda0 = MIN(zda0, zdamax) 364 268 zetamax = MIN( zdh0, hR(ii,ij,klbnd) ) - hL(ii,ij,klbnd) 269 270 IF( zetamax > 0.0 ) THEN 271 zx1 = zetamax 272 zx2 = 0.5 * zetamax * zetamax 273 zda0 = g1(ii,ij,klbnd) * zx2 + g0(ii,ij,klbnd) * zx1 ! ice area removed 274 zdamax = a_i(ii,ij,klbnd) * (1.0 - ht_i(ii,ij,klbnd) / zht_i_b(ii,ij,klbnd) ) ! Constrain new thickness <= ht_i 275 zda0 = MIN( zda0, zdamax ) ! ice area lost due to melting 276 ! of thin ice (zdamax > 0) 365 277 ! Remove area, conserving volume 366 ht_i(ii,ij,klbnd) = ht_i(ii,ij,klbnd) & 367 * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 ) 278 ht_i(ii,ij,klbnd) = ht_i(ii,ij,klbnd) * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 ) 368 279 a_i(ii,ij,klbnd) = a_i(ii,ij,klbnd) - zda0 369 v_i(ii,ij,klbnd) = a_i(ii,ij,klbnd)*ht_i(ii,ij,klbnd) ! clem-useless ? 370 ENDIF ! zetamax > 0 371 ! ji, a_i > epsi10 372 373 ELSE ! if ice accretion 374 ! ji, a_i > epsi10; zdh0 > 0 375 zhbnew(ii,ij,klbnd-1) = MIN(zdh0,hi_max(klbnd)) 376 ! zhbnew was 0, and is shifted to the right to account for thin ice 377 ! growth in openwater (F0 = f1) 378 ENDIF ! zdh0 379 380 ! a_i > epsi10 381 ENDIF ! a_i > epsi10 382 383 END DO ! ji 280 v_i(ii,ij,klbnd) = a_i(ii,ij,klbnd) * ht_i(ii,ij,klbnd) ! clem-useless ? 281 ENDIF 282 283 ELSE ! if ice accretion zdh0 > 0 284 ! zhbnew was 0, and is shifted to the right to account for thin ice growth in openwater (F0 = f1) 285 zhbnew(ii,ij,klbnd-1) = MIN( zdh0, hi_max(klbnd) ) 286 ENDIF 287 288 ENDIF 289 290 END DO 384 291 385 292 !- 7.3 g(h) for each thickness category 386 293 DO jl = klbnd, kubnd 387 CALL lim_itd_fitline( jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl),&388 g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag)294 CALL lim_itd_fitline( jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl), & 295 & g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag ) 389 296 END DO 390 297 … … 406 313 ij = nind_j(ji) 407 314 408 IF (zhbnew(ii,ij,jl) .gt. hi_max(jl)) THEN ! transfer from jl to jl+1 409 315 IF (zhbnew(ii,ij,jl) > hi_max(jl)) THEN ! transfer from jl to jl+1 410 316 ! left and right integration limits in eta space 411 zvetamin(ji) = MAX( hi_max(jl), hL(ii,ij,jl)) - hL(ii,ij,jl)412 zvetamax(ji) = MIN( zhbnew(ii,ij,jl), hR(ii,ij,jl)) - hL(ii,ij,jl)317 zvetamin(ji) = MAX( hi_max(jl), hL(ii,ij,jl) ) - hL(ii,ij,jl) 318 zvetamax(ji) = MIN( zhbnew(ii,ij,jl), hR(ii,ij,jl) ) - hL(ii,ij,jl) 413 319 zdonor(ii,ij,jl) = jl 414 320 415 ELSE ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl 416 321 ELSE ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl 417 322 ! left and right integration limits in eta space 418 323 zvetamin(ji) = 0.0 419 zvetamax(ji) = MIN( hi_max(jl), hR(ii,ij,jl+1)) - hL(ii,ij,jl+1)324 zvetamax(ji) = MIN( hi_max(jl), hR(ii,ij,jl+1) ) - hL(ii,ij,jl+1) 420 325 zdonor(ii,ij,jl) = jl + 1 421 326 422 ENDIF ! zhbnew(jl) > hi_max(jl)423 424 zetamax = MAX( zvetamax(ji), zvetamin(ji)) ! no transfer if etamax < etamin327 ENDIF 328 329 zetamax = MAX( zvetamax(ji), zvetamin(ji) ) ! no transfer if etamax < etamin 425 330 zetamin = zvetamin(ji) 426 331 427 332 zx1 = zetamax - zetamin 428 zwk1 = zetamin *zetamin429 zwk2 = zetamax *zetamax430 zx2 = 0.5 * ( zwk2 - zwk1)333 zwk1 = zetamin * zetamin 334 zwk2 = zetamax * zetamax 335 zx2 = 0.5 * ( zwk2 - zwk1 ) 431 336 zwk1 = zwk1 * zetamin 432 337 zwk2 = zwk2 * zetamax 433 zx3 = 1.0 /3.0 * (zwk2 - zwk1)338 zx3 = 1.0 / 3.0 * ( zwk2 - zwk1 ) 434 339 nd = zdonor(ii,ij,jl) 435 340 zdaice(ii,ij,jl) = g1(ii,ij,nd)*zx2 + g0(ii,ij,nd)*zx1 436 341 zdvice(ii,ij,jl) = g1(ii,ij,nd)*zx3 + g0(ii,ij,nd)*zx2 + zdaice(ii,ij,jl)*hL(ii,ij,nd) 437 342 438 END DO ! ji439 END DO ! jl klbnd -> kubnd - 1343 END DO 344 END DO 440 345 441 346 !!---------------------------------------------------------------------------------------------- … … 451 356 ii = nind_i(ji) 452 357 ij = nind_j(ji) 453 IF ( a_i(ii,ij,1) > epsi10 .AND. ht_i(ii,ij,1) < hiclim) THEN454 a_i (ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / hiclim455 ht_i(ii,ij,1) = hiclim358 IF ( a_i(ii,ij,1) > epsi10 .AND. ht_i(ii,ij,1) < rn_himin ) THEN 359 a_i (ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / rn_himin 360 ht_i(ii,ij,1) = rn_himin 456 361 ENDIF 457 END DO !ji362 END DO 458 363 459 364 !!---------------------------------------------------------------------------------------------- … … 479 384 ENDIF 480 385 481 CALL wrk_dealloc( jpi,jpj, zremap_flag ) ! integer482 CALL wrk_dealloc( jpi,jpj,jpl-1, zdonor ) ! integer386 CALL wrk_dealloc( jpi,jpj, zremap_flag ) 387 CALL wrk_dealloc( jpi,jpj,jpl-1, zdonor ) 483 388 CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es ) 484 389 CALL wrk_dealloc( jpi,jpj,jpl-1, zdaice, zdvice ) 485 390 CALL wrk_dealloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 ) 486 391 CALL wrk_dealloc( (jpi+1)*(jpj+1), zvetamin, zvetamax ) 487 CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) ! integer392 CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 488 393 CALL wrk_dealloc( jpi,jpj, zhb0,zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final ) 489 394 … … 491 396 492 397 493 SUBROUTINE lim_itd_fitline( num_cat, HbL, Hbr, hice, & 494 & g0, g1, hL, hR, zremap_flag ) 398 SUBROUTINE lim_itd_fitline( num_cat, HbL, Hbr, hice, g0, g1, hL, hR, zremap_flag ) 495 399 !!------------------------------------------------------------------ 496 400 !! *** ROUTINE lim_itd_fitline *** … … 511 415 INTEGER , DIMENSION(jpi,jpj), INTENT(in ) :: zremap_flag ! 512 416 ! 513 INTEGER :: ji,jj! horizontal indices417 INTEGER :: ji,jj ! horizontal indices 514 418 REAL(wp) :: zh13 ! HbL + 1/3 * (HbR - HbL) 515 419 REAL(wp) :: zh23 ! HbL + 2/3 * (HbR - HbL) … … 518 422 !!------------------------------------------------------------------ 519 423 ! 520 !521 424 DO jj = 1, jpj 522 425 DO ji = 1, jpi 523 426 ! 524 427 IF( zremap_flag(ji,jj) == 1 .AND. a_i(ji,jj,num_cat) > epsi10 & 525 & .AND. hice(ji,jj) > 0._wp )THEN428 & .AND. hice(ji,jj) > 0._wp ) THEN 526 429 527 430 ! Initialize hL and hR 528 529 431 hL(ji,jj) = HbL(ji,jj) 530 432 hR(ji,jj) = HbR(ji,jj) 531 433 532 434 ! Change hL or hR if hice falls outside central third of range 533 534 zh13 = 1.0/3.0 * (2.0*hL(ji,jj) + hR(ji,jj)) 535 zh23 = 1.0/3.0 * (hL(ji,jj) + 2.0*hR(ji,jj)) 435 zh13 = 1.0 / 3.0 * ( 2.0 * hL(ji,jj) + hR(ji,jj) ) 436 zh23 = 1.0 / 3.0 * ( hL(ji,jj) + 2.0 * hR(ji,jj) ) 536 437 537 438 IF ( hice(ji,jj) < zh13 ) THEN ; hR(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hL(ji,jj) … … 540 441 541 442 ! Compute coefficients of g(eta) = g0 + g1*eta 542 543 443 zdhr = 1._wp / (hR(ji,jj) - hL(ji,jj)) 544 444 zwk1 = 6._wp * a_i(ji,jj,num_cat) * zdhr 545 445 zwk2 = ( hice(ji,jj) - hL(ji,jj) ) * zdhr 546 g0(ji,jj) = zwk1 * ( 2._wp /3._wp - zwk2 )547 g1(ji,jj) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5)446 g0(ji,jj) = zwk1 * ( 2._wp / 3._wp - zwk2 ) 447 g1(ji,jj) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5 ) 548 448 ! 549 ELSE 449 ELSE ! remap_flag = .false. or a_i < epsi10 550 450 hL(ji,jj) = 0._wp 551 451 hR(ji,jj) = 0._wp 552 452 g0(ji,jj) = 0._wp 553 453 g1(ji,jj) = 0._wp 554 ENDIF ! a_i > epsi10454 ENDIF 555 455 ! 556 456 END DO … … 576 476 577 477 INTEGER :: ji, jj, jl, jl2, jl1, jk ! dummy loop indices 578 INTEGER :: ii, ij ! indices when changing from 2D-1D is done478 INTEGER :: ii, ij ! indices when changing from 2D-1D is done 579 479 580 480 REAL(wp), POINTER, DIMENSION(:,:,:) :: zaTsfn … … 589 489 INTEGER, POINTER, DIMENSION(:) :: nind_i, nind_j ! compressed indices for i/j directions 590 490 591 INTEGER :: nbrem ! number of cells with ice to transfer 592 593 LOGICAL :: zdaice_negative ! true if daice < -puny 594 LOGICAL :: zdvice_negative ! true if dvice < -puny 595 LOGICAL :: zdaice_greater_aicen ! true if daice > aicen 596 LOGICAL :: zdvice_greater_vicen ! true if dvice > vicen 491 INTEGER :: nbrem ! number of cells with ice to transfer 597 492 !!------------------------------------------------------------------ 598 493 599 494 CALL wrk_alloc( jpi,jpj,jpl, zaTsfn ) 600 495 CALL wrk_alloc( jpi,jpj, zworka ) 601 CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) ! integer496 CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 602 497 603 498 !---------------------------------------------------------------------------------------------- … … 606 501 607 502 DO jl = klbnd, kubnd 608 zaTsfn(:,:,jl) = a_i(:,:,jl)*t_su(:,:,jl) 609 END DO 610 611 !---------------------------------------------------------------------------------------------- 612 ! 2) Check for daice or dvice out of range, allowing for roundoff error 613 !---------------------------------------------------------------------------------------------- 614 ! Note: zdaice < 0 or zdvice < 0 usually happens when category jl 615 ! has a small area, with h(n) very close to a boundary. Then 616 ! the coefficients of g(h) are large, and the computed daice and 617 ! dvice can be in error. If this happens, it is best to transfer 618 ! either the entire category or nothing at all, depending on which 619 ! side of the boundary hice(n) lies. 620 !----------------------------------------------------------------- 621 DO jl = klbnd, kubnd-1 622 623 zdaice_negative = .false. 624 zdvice_negative = .false. 625 zdaice_greater_aicen = .false. 626 zdvice_greater_vicen = .false. 627 628 DO jj = 1, jpj 629 DO ji = 1, jpi 630 631 IF (zdonor(ji,jj,jl) .GT. 0) THEN 632 jl1 = zdonor(ji,jj,jl) 633 634 IF (zdaice(ji,jj,jl) .LT. 0.0) THEN 635 IF (zdaice(ji,jj,jl) .GT. -epsi10) THEN 636 IF ( ( jl1.EQ.jl .AND. ht_i(ji,jj,jl1) .GT. hi_max(jl) ) & 637 .OR. & 638 ( jl1.EQ.jl+1 .AND. ht_i(ji,jj,jl1) .LE. hi_max(jl) ) & 639 ) THEN 640 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category 641 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 642 ELSE 643 zdaice(ji,jj,jl) = 0.0 ! shift no ice 644 zdvice(ji,jj,jl) = 0.0 645 ENDIF 646 ELSE 647 zdaice_negative = .true. 648 ENDIF 649 ENDIF 650 651 IF (zdvice(ji,jj,jl) .LT. 0.0) THEN 652 IF (zdvice(ji,jj,jl) .GT. -epsi10 ) THEN 653 IF ( ( jl1.EQ.jl .AND. ht_i(ji,jj,jl1).GT.hi_max(jl) ) & 654 .OR. & 655 ( jl1.EQ.jl+1 .AND. ht_i(ji,jj,jl1) .LE. hi_max(jl) ) & 656 ) THEN 657 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category 658 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 659 ELSE 660 zdaice(ji,jj,jl) = 0.0 ! shift no ice 661 zdvice(ji,jj,jl) = 0.0 662 ENDIF 663 ELSE 664 zdvice_negative = .true. 665 ENDIF 666 ENDIF 667 668 ! If daice is close to aicen, set daice = aicen. 669 IF (zdaice(ji,jj,jl) .GT. a_i(ji,jj,jl1) - epsi10 ) THEN 670 IF (zdaice(ji,jj,jl) .LT. a_i(ji,jj,jl1)+epsi10) THEN 671 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 672 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 673 ELSE 674 zdaice_greater_aicen = .true. 675 ENDIF 676 ENDIF 677 678 IF (zdvice(ji,jj,jl) .GT. v_i(ji,jj,jl1)-epsi10) THEN 679 IF (zdvice(ji,jj,jl) .LT. v_i(ji,jj,jl1)+epsi10) THEN 680 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 681 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 682 ELSE 683 zdvice_greater_vicen = .true. 684 ENDIF 685 ENDIF 686 687 ENDIF ! donor > 0 688 END DO ! i 689 END DO ! j 690 691 END DO !jl 503 zaTsfn(:,:,jl) = a_i(:,:,jl) * t_su(:,:,jl) 504 END DO 692 505 693 506 !------------------------------------------------------------------------------- 694 ! 3) Transfer volume and energy between categories507 ! 2) Transfer volume and energy between categories 695 508 !------------------------------------------------------------------------------- 696 509 … … 699 512 DO jj = 1, jpj 700 513 DO ji = 1, jpi 701 IF (zdaice(ji,jj,jl) .GT.0.0 ) THEN ! daice(n) can be < puny514 IF (zdaice(ji,jj,jl) > 0.0 ) THEN ! daice(n) can be < puny 702 515 nbrem = nbrem + 1 703 516 nind_i(nbrem) = ji 704 517 nind_j(nbrem) = jj 705 ENDIF ! tmask518 ENDIF 706 519 END DO 707 520 END DO … … 712 525 713 526 jl1 = zdonor(ii,ij,jl) 714 rswitch = MAX( 0.0 , SIGN( 1.0, v_i(ii,ij,jl1) - epsi10 ) )715 zworka(ii,ij) = zdvice(ii,ij,jl) / MAX(v_i(ii,ij,jl1),epsi10) * rswitch527 rswitch = MAX( 0._wp , SIGN( 1._wp , v_i(ii,ij,jl1) - epsi10 ) ) 528 zworka(ii,ij) = zdvice(ii,ij,jl) / MAX( v_i(ii,ij,jl1), epsi10 ) * rswitch 716 529 IF( jl1 == jl) THEN ; jl2 = jl1+1 717 ELSE 530 ELSE ; jl2 = jl 718 531 ENDIF 719 532 … … 721 534 ! Ice areas 722 535 !-------------- 723 724 536 a_i(ii,ij,jl1) = a_i(ii,ij,jl1) - zdaice(ii,ij,jl) 725 537 a_i(ii,ij,jl2) = a_i(ii,ij,jl2) + zdaice(ii,ij,jl) … … 728 540 ! Ice volumes 729 541 !-------------- 730 731 542 v_i(ii,ij,jl1) = v_i(ii,ij,jl1) - zdvice(ii,ij,jl) 732 543 v_i(ii,ij,jl2) = v_i(ii,ij,jl2) + zdvice(ii,ij,jl) … … 735 546 ! Snow volumes 736 547 !-------------- 737 738 548 zdvsnow = v_s(ii,ij,jl1) * zworka(ii,ij) 739 549 v_s(ii,ij,jl1) = v_s(ii,ij,jl1) - zdvsnow … … 743 553 ! Snow heat content 744 554 !-------------------- 745 746 555 zdesnow = e_s(ii,ij,1,jl1) * zworka(ii,ij) 747 556 e_s(ii,ij,1,jl1) = e_s(ii,ij,1,jl1) - zdesnow … … 751 560 ! Ice age 752 561 !-------------- 753 754 562 zdo_aice = oa_i(ii,ij,jl1) * zdaice(ii,ij,jl) 755 563 oa_i(ii,ij,jl1) = oa_i(ii,ij,jl1) - zdo_aice … … 759 567 ! Ice salinity 760 568 !-------------- 761 762 569 zdsm_vice = smv_i(ii,ij,jl1) * zworka(ii,ij) 763 570 smv_i(ii,ij,jl1) = smv_i(ii,ij,jl1) - zdsm_vice … … 767 574 ! Surface temperature 768 575 !--------------------- 769 770 576 zdaTsf = t_su(ii,ij,jl1) * zdaice(ii,ij,jl) 771 577 zaTsfn(ii,ij,jl1) = zaTsfn(ii,ij,jl1) - zdaTsf 772 578 zaTsfn(ii,ij,jl2) = zaTsfn(ii,ij,jl2) + zdaTsf 773 579 774 END DO ! ji580 END DO 775 581 776 582 !------------------ … … 779 585 780 586 DO jk = 1, nlay_i 781 !CDIR NODEP782 587 DO ji = 1, nbrem 783 588 ii = nind_i(ji) … … 785 590 786 591 jl1 = zdonor(ii,ij,jl) 787 IF (jl1 .EQ.jl) THEN592 IF (jl1 == jl) THEN 788 593 jl2 = jl+1 789 594 ELSE ! n1 = n+1 … … 794 599 e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - zdeice 795 600 e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + zdeice 796 END DO ! ji797 END DO ! jk601 END DO 602 END DO 798 603 799 604 END DO ! boundaries, 1 to ncat-1 … … 809 614 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl) 810 615 t_su(ji,jj,jl) = zaTsfn(ji,jj,jl) / a_i(ji,jj,jl) 811 rswitch = 1.0 - MAX(0.0,SIGN(1.0,-v_s(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes812 616 ELSE 813 617 ht_i(ji,jj,jl) = 0._wp 814 t_su(ji,jj,jl) = rt t618 t_su(ji,jj,jl) = rt0 815 619 ENDIF 816 END DO ! ji817 END DO ! jj818 END DO ! jl620 END DO 621 END DO 622 END DO 819 623 ! 820 624 CALL wrk_dealloc( jpi,jpj,jpl, zaTsfn ) 821 625 CALL wrk_dealloc( jpi,jpj, zworka ) 822 CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) ! integer626 CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 823 627 ! 824 628 END SUBROUTINE lim_itd_shiftice … … 846 650 REAL(wp), POINTER, DIMENSION(:,:) :: vt_s_init, vt_s_final ! snow volume summed over categories 847 651 !!------------------------------------------------------------------ 848 !! clem 2014/04: be carefull, rebining does not conserve salt(maybe?) => the difference is taken into account in limupdate849 652 850 653 CALL wrk_alloc( jpi,jpj,jpl, zdonor ) ! interger … … 864 667 DO jj = 1, jpj 865 668 DO ji = 1, jpi 866 IF( a_i(ji,jj,jl) > epsi10 ) THEN 867 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 868 ELSE 869 ht_i(ji,jj,jl) = 0._wp 870 ENDIF 669 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) ) 670 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch 871 671 END DO 872 672 END DO … … 874 674 875 675 !------------------------------------------------------------------------------ 876 ! 2) Make sure thickness of cat klbnd is at least hi_max(klbnd) 877 !------------------------------------------------------------------------------ 878 DO jj = 1, jpj 879 DO ji = 1, jpi 880 IF( a_i(ji,jj,klbnd) > epsi10 ) THEN 881 IF( ht_i(ji,jj,klbnd) <= hi_max(0) .AND. hi_max(0) > 0._wp ) THEN 882 a_i(ji,jj,klbnd) = v_i(ji,jj,klbnd) / hi_max(0) 883 ht_i(ji,jj,klbnd) = hi_max(0) 884 ENDIF 885 ENDIF 886 END DO 887 END DO 888 889 !------------------------------------------------------------------------------ 890 ! 3) If a category thickness is not in bounds, shift the 676 ! 2) If a category thickness is not in bounds, shift the 891 677 ! entire area, volume, and energy to the neighboring category 892 678 !------------------------------------------------------------------------------ … … 917 703 zdonor(ji,jj,jl) = jl 918 704 ! begin TECLIM change 919 !zdaice(ji,jj,jl) = a_i(ji,jj,jl)920 !zdvice(ji,jj,jl) = v_i(ji,jj,jl)921 705 !zdaice(ji,jj,jl) = a_i(ji,jj,jl) * 0.5_wp 922 706 !zdvice(ji,jj,jl) = v_i(ji,jj,jl)-zdaice(ji,jj,jl)*(hi_max(jl)+hi_max(jl-1)) * 0.5_wp 923 707 ! end TECLIM change 924 708 ! clem: how much of a_i you send in cat sup is somewhat arbitrary 925 zdaice(ji,jj,jl) = a_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - hi_max(jl) + epsi 10 ) / ht_i(ji,jj,jl)926 zdvice(ji,jj,jl) = v_i(ji,jj,jl) - ( a_i(ji,jj,jl) - zdaice(ji,jj,jl) ) * ( hi_max(jl) - epsi 10 )709 zdaice(ji,jj,jl) = a_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - hi_max(jl) + epsi20 ) / ht_i(ji,jj,jl) 710 zdvice(ji,jj,jl) = v_i(ji,jj,jl) - ( a_i(ji,jj,jl) - zdaice(ji,jj,jl) ) * ( hi_max(jl) - epsi20 ) 927 711 ENDIF 928 END DO ! ji929 END DO ! jj712 END DO 713 END DO 930 714 IF(lk_mpp) CALL mpp_max( zshiftflag ) 931 715 … … 938 722 ENDIF 939 723 ! 940 END DO ! jl724 END DO 941 725 942 726 !---------------------------- … … 951 735 zshiftflag = 0 952 736 953 !clem-change954 737 DO jj = 1, jpj 955 738 DO ji = 1, jpi … … 961 744 zdvice(ji,jj,jl) = v_i(ji,jj,jl+1) 962 745 ENDIF 963 END DO ! ji964 END DO ! jj746 END DO 747 END DO 965 748 966 749 IF(lk_mpp) CALL mpp_max( zshiftflag ) … … 973 756 zdvice(:,:,jl) = 0._wp 974 757 ENDIF 975 !clem-change 976 977 ! ! clem-change begin: why not doing that? 978 ! DO jj = 1, jpj 979 ! DO ji = 1, jpi 980 ! IF( a_i(ji,jj,jl+1) > epsi10 .AND. ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN 981 ! ht_i(ji,jj,jl+1) = hi_max(jl) + epsi10 982 ! a_i (ji,jj,jl+1) = v_i(ji,jj,jl+1) / ht_i(ji,jj,jl+1) 983 ! ENDIF 984 ! END DO ! ji 985 ! END DO ! jj 986 ! clem-change end 987 988 END DO ! jl 758 759 END DO 989 760 990 761 !------------------------------------------------------------------------------ 991 ! 4) Conservation check762 ! 3) Conservation check 992 763 !------------------------------------------------------------------------------ 993 764 … … 1002 773 ENDIF 1003 774 ! 1004 CALL wrk_dealloc( jpi,jpj,jpl, zdonor ) ! interger775 CALL wrk_dealloc( jpi,jpj,jpl, zdonor ) 1005 776 CALL wrk_dealloc( jpi,jpj,jpl, zdaice, zdvice ) 1006 777 CALL wrk_dealloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final ) … … 1013 784 !!---------------------------------------------------------------------- 1014 785 CONTAINS 1015 SUBROUTINE lim_itd_th ! Empty routines1016 END SUBROUTINE lim_itd_th1017 SUBROUTINE lim_itd_th_ini1018 END SUBROUTINE lim_itd_th_ini1019 786 SUBROUTINE lim_itd_th_rem 1020 787 END SUBROUTINE lim_itd_th_rem
Note: See TracChangeset
for help on using the changeset viewer.
