Changeset 5034 for branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
- Timestamp:
- 2015-01-15T14:48:42+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
r4333 r5034 6 6 !! History : - ! (W. H. Lipscomb and E.C. Hunke) CICE (c) original code 7 7 !! 3.0 ! 2005-12 (M. Vancoppenolle) adaptation to LIM-3 8 !! - ! 2006-06 (M. Vancoppenolle) adaptation to include salt, age and types8 !! - ! 2006-06 (M. Vancoppenolle) adaptation to include salt, age 9 9 !! - ! 2007-04 (M. Vancoppenolle) Mass conservation checked 10 10 !!---------------------------------------------------------------------- … … 35 35 USE lib_fortran ! to use key_nosignedzero 36 36 USE timing ! Timing 37 USE limcons ! conservation tests 37 38 38 39 IMPLICIT NONE … … 44 45 PUBLIC lim_itd_fitline 45 46 PUBLIC lim_itd_shiftice 46 47 REAL(wp) :: epsi10 = 1.e-10_wp !48 REAL(wp) :: epsi06 = 1.e-6_wp !49 47 50 48 !!---------------------------------------------------------------------- … … 65 63 INTEGER, INTENT(in) :: kt ! time step index 66 64 ! 67 INTEGER :: j l, ja, jm, jbnd1, jbnd2 ! ice typesdummy loop index68 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)69 REAL(wp) :: z chk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset)65 INTEGER :: ji, jj, jk, jl ! dummy loop index 66 ! 67 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 70 68 !!------------------------------------------------------------------ 71 69 IF( nn_timing == 1 ) CALL timing_start('limitd_th') 72 70 73 ! ------------------------------- 74 !- check conservation (C Rousset) 75 IF (ln_limdiahsb) THEN 76 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 77 zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 78 zchk_fw_b = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 79 zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 80 ENDIF 81 !- check conservation (C Rousset) 82 ! ------------------------------- 71 ! conservation test 72 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_th', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 83 73 84 74 IF( kt == nit000 .AND. lwp ) THEN … … 93 83 ! Given thermodynamic growth rates, transport ice between 94 84 ! thickness categories. 95 DO jm = 1, jpm 96 jbnd1 = ice_cat_bounds(jm,1) 97 jbnd2 = ice_cat_bounds(jm,2) 98 IF( ice_ncat_types(jm) > 1 ) CALL lim_itd_th_rem( jbnd1, jbnd2, jm, kt ) 99 END DO 85 IF( jpl > 1 ) CALL lim_itd_th_rem( 1, jpl, kt ) 100 86 ! 101 87 CALL lim_var_glo2eqv ! only for info … … 105 91 ! 3) Add frazil ice growing in leads. 106 92 !------------------------------------------------------------------------------| 107 108 93 CALL lim_thd_lac 109 94 CALL lim_var_glo2eqv ! only for info 110 111 IF(ln_ctl) THEN ! Control print95 96 IF(ln_ctl) THEN ! Control print 112 97 CALL prt_ctl_info(' ') 113 98 CALL prt_ctl_info(' - Cell values : ') … … 131 116 CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_itd_th : sm_i : ') 132 117 CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_itd_th : smv_i : ') 133 DO j a= 1, nlay_i118 DO jk = 1, nlay_i 134 119 CALL prt_ctl_info(' ') 135 CALL prt_ctl_info(' - Layer : ', ivar1=j a)120 CALL prt_ctl_info(' - Layer : ', ivar1=jk) 136 121 CALL prt_ctl_info(' ~~~~~~~') 137 CALL prt_ctl(tab2d_1=t_i(:,:,j a,jl) , clinfo1= ' lim_itd_th : t_i : ')138 CALL prt_ctl(tab2d_1=e_i(:,:,j a,jl) , clinfo1= ' lim_itd_th : e_i : ')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 : ') 139 124 END DO 140 125 END DO 141 126 ENDIF 142 127 ! 143 ! ------------------------------- 144 !- check conservation (C Rousset) 145 IF( ln_limdiahsb ) THEN 146 zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 147 zchk_fw = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 148 149 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice 150 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 151 152 zchk_vmin = glob_min(v_i) 153 zchk_amax = glob_max(SUM(a_i,dim=3)) 154 zchk_amin = glob_min(a_i) 155 156 IF(lwp) THEN 157 IF ( ABS( zchk_v_i ) > 1.e-5 ) WRITE(numout,*) 'violation volume [m3/day] (limitd_th) = ',(zchk_v_i * rday) 158 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_th) = ',(zchk_smv * rday) 159 IF ( zchk_vmin < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limitd_th) = ',(zchk_vmin * 1.e-3) 160 IF ( zchk_amax > amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax (limitd_th) = ',zchk_amax 161 IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limitd_th) = ',zchk_amin 162 ENDIF 163 ENDIF 164 !- check conservation (C Rousset) 165 ! ------------------------------- 128 ! conservation test 129 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_th', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 166 130 ! 167 131 IF( nn_timing == 1 ) CALL timing_stop('limitd_th') … … 169 133 ! 170 134 171 SUBROUTINE lim_itd_th_rem( klbnd, kubnd, ntyp,kt )135 SUBROUTINE lim_itd_th_rem( klbnd, kubnd, kt ) 172 136 !!------------------------------------------------------------------ 173 137 !! *** ROUTINE lim_itd_th_rem *** … … 182 146 INTEGER , INTENT (in) :: klbnd ! Start thickness category index point 183 147 INTEGER , INTENT (in) :: kubnd ! End point on which the the computation is applied 184 INTEGER , INTENT (in) :: ntyp ! Number of the type used185 148 INTEGER , INTENT (in) :: kt ! Ocean time step 186 149 ! … … 190 153 REAL(wp) :: zx1, zwk1, zdh0, zetamin, zdamax ! local scalars 191 154 REAL(wp) :: zx2, zwk2, zda0, zetamax ! - - 192 REAL(wp) :: zx3, zareamin , zindb! - -155 REAL(wp) :: zx3, zareamin ! - - 193 156 CHARACTER (len = 15) :: fieldid 194 157 … … 200 163 REAL(wp), POINTER, DIMENSION(:,:,:) :: hL ! left boundary for the ITD for each thickness 201 164 REAL(wp), POINTER, DIMENSION(:,:,:) :: hR ! left boundary for the ITD for each thickness 202 REAL(wp), POINTER, DIMENSION(:,:,:) :: zht_i_ o! old ice thickness165 REAL(wp), POINTER, DIMENSION(:,:,:) :: zht_i_b ! old ice thickness 203 166 REAL(wp), POINTER, DIMENSION(:,:,:) :: dummy_es 204 167 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdaice, zdvice ! local increment of ice area and volume … … 218 181 CALL wrk_alloc( jpi,jpj, zremap_flag ) ! integer 219 182 CALL wrk_alloc( jpi,jpj,jpl-1, zdonor ) ! integer 220 CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_ o, dummy_es )183 CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es ) 221 184 CALL wrk_alloc( jpi,jpj,jpl-1, zdaice, zdvice ) 222 185 CALL wrk_alloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 ) … … 247 210 WRITE(numout,*) ' klbnd : ', klbnd 248 211 WRITE(numout,*) ' kubnd : ', kubnd 249 WRITE(numout,*) ' ntyp : ', ntyp250 212 ENDIF 251 213 … … 254 216 DO jj = 1, jpj 255 217 DO ji = 1, jpi 256 zindb= 1.0 - MAX( 0.0, SIGN( 1.0, - a_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes257 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zindb258 zindb = 1.0 - MAX( 0.0, SIGN( 1.0, - old_a_i(ji,jj,jl) + epsi10) ) !0 if no ice and 1 if yes259 zht_i_ o(ji,jj,jl) = old_v_i(ji,jj,jl) / MAX( old_a_i(ji,jj,jl), epsi10 ) * zindb260 IF( a_i(ji,jj,jl) > epsi 06 ) zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_o(ji,jj,jl)218 rswitch = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes 219 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 yes 221 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) 261 223 END DO 262 224 END DO … … 302 264 ij = nind_j(ji) 303 265 ! 304 IF ( ( zht_i_o(ii,ij,jl) .GT. epsi10 ) .AND. &305 ( zht_i_o(ii,ij,jl+1) .GT. epsi10 )) THEN266 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 ) THEN 306 268 !interpolate between adjacent category growth rates 307 zslope = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / & 308 ( zht_i_o (ii,ij,jl+1) - zht_i_o (ii,ij,jl) ) 309 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + & 310 zslope * ( hi_max(jl) - zht_i_o(ii,ij,jl) ) 311 ELSEIF (zht_i_o(ii,ij,jl).gt.epsi10) THEN 269 zslope = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_b(ii,ij,jl+1) - zht_i_b(ii,ij,jl) ) 270 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_b(ii,ij,jl) ) 271 ELSEIF ( a_i_b(ii,ij,jl) > epsi10) THEN 312 272 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) 313 ELSEIF ( zht_i_o(ii,ij,jl+1).gt.epsi10) THEN273 ELSEIF ( a_i_b(ii,ij,jl+1) > epsi10) THEN 314 274 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1) 315 ELSE316 zhbnew(ii,ij,jl) = hi_max(jl)317 275 ENDIF 318 276 END DO … … 320 278 !- 4.2 Check that each zhbnew lies between adjacent values of ice thickness 321 279 DO ji = 1, nbrem 322 ! jl, ji323 280 ii = nind_i(ji) 324 281 ij = nind_j(ji) 325 ! jl, ji 326 IF ( ( a_i(ii,ij,jl) .GT.epsi10) .AND. & 327 ( ht_i(ii,ij,jl).GE. zhbnew(ii,ij,jl) ) & 328 ) THEN 282 IF( a_i(ii,ij,jl) > epsi10 .AND. ht_i(ii,ij,jl) >= zhbnew(ii,ij,jl) ) THEN 329 283 zremap_flag(ii,ij) = 0 330 ELSEIF ( ( a_i(ii,ij,jl+1) .GT. epsi10 ) .AND. & 331 ( ht_i(ii,ij,jl+1).LE. zhbnew(ii,ij,jl) ) & 332 ) THEN 284 ELSEIF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) <= zhbnew(ii,ij,jl) ) THEN 333 285 zremap_flag(ii,ij) = 0 334 286 ENDIF 335 287 336 288 !- 4.3 Check that each zhbnew does not exceed maximal values hi_max 337 ! jl, ji 338 IF (zhbnew(ii,ij,jl).gt.hi_max(jl+1)) THEN 339 zremap_flag(ii,ij) = 0 340 ENDIF 341 ! jl, ji 342 IF (zhbnew(ii,ij,jl).lt.hi_max(jl-1)) THEN 343 zremap_flag(ii,ij) = 0 344 ENDIF 345 ! jl, ji 346 END DO !ji 347 ! ji 289 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 348 293 END DO !jl 349 294 … … 354 299 DO jj = 1, jpj 355 300 DO ji = 1, jpi 356 IF 301 IF( zremap_flag(ji,jj) == 1 ) THEN 357 302 nbrem = nbrem + 1 358 303 nind_i(nbrem) = ji 359 304 nind_j(nbrem) = jj 360 305 ENDIF 361 END DO !ji362 END DO !jj306 END DO 307 END DO 363 308 364 309 !----------------------------------------------------------------------------------------------- … … 367 312 DO jj = 1, jpj 368 313 DO ji = 1, jpi 369 zhb0(ji,jj) = hi_max _typ(0,ntyp) ! 0eme370 zhb1(ji,jj) = hi_max _typ(1,ntyp) ! 1er314 zhb0(ji,jj) = hi_max(0) ! 0eme 315 zhb1(ji,jj) = hi_max(1) ! 1er 371 316 372 317 zhbnew(ji,jj,klbnd-1) = 0._wp … … 380 325 ENDIF 381 326 382 IF( zhbnew(ji,jj,kubnd) < hi_max(kubnd-1) ) 327 IF( zhbnew(ji,jj,kubnd) < hi_max(kubnd-1) ) zhbnew(ji,jj,kubnd) = hi_max(kubnd-1) 383 328 384 329 END DO !jj … … 389 334 !----------------------------------------------------------------------------------------------- 390 335 !- 7.1 g(h) for category 1 at start of time step 391 CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_ o(:,:,klbnd), &336 CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd), & 392 337 & g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd), & 393 338 & hR(:,:,klbnd), zremap_flag ) … … 414 359 ! Constrain new thickness <= ht_i 415 360 zdamax = a_i(ii,ij,klbnd) * & 416 (1.0 - ht_i(ii,ij,klbnd)/zht_i_ o(ii,ij,klbnd)) ! zdamax > 0361 (1.0 - ht_i(ii,ij,klbnd)/zht_i_b(ii,ij,klbnd)) ! zdamax > 0 417 362 !ice area lost due to melting of thin ice 418 363 zda0 = MIN(zda0, zdamax) … … 428 373 ELSE ! if ice accretion 429 374 ! ji, a_i > epsi10; zdh0 > 0 430 IF ( ntyp .EQ. 1 )zhbnew(ii,ij,klbnd-1) = MIN(zdh0,hi_max(klbnd))375 zhbnew(ii,ij,klbnd-1) = MIN(zdh0,hi_max(klbnd)) 431 376 ! zhbnew was 0, and is shifted to the right to account for thin ice 432 377 ! growth in openwater (F0 = f1) 433 IF ( ntyp .NE. 1 ) zhbnew(ii,ij,0) = 0434 ! in other types there is435 ! no open water growth (F0 = 0)436 378 ENDIF ! zdh0 437 379 … … 444 386 DO jl = klbnd, kubnd 445 387 CALL lim_itd_fitline(jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl), & 446 g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), & 447 zremap_flag) 388 g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag) 448 389 END DO 449 390 … … 493 434 nd = zdonor(ii,ij,jl) 494 435 zdaice(ii,ij,jl) = g1(ii,ij,nd)*zx2 + g0(ii,ij,nd)*zx1 495 zdvice(ii,ij,jl) = g1(ii,ij,nd)*zx3 + g0(ii,ij,nd)*zx2 + & 496 zdaice(ii,ij,jl)*hL(ii,ij,nd) 436 zdvice(ii,ij,jl) = g1(ii,ij,nd)*zx3 + g0(ii,ij,nd)*zx2 + zdaice(ii,ij,jl)*hL(ii,ij,nd) 497 437 498 438 END DO ! ji … … 511 451 ii = nind_i(ji) 512 452 ij = nind_j(ji) 513 IF ( ( a_i(ii,ij,1) > epsi10 ) .AND. ( ht_i(ii,ij,1) < hiclim )) THEN453 IF ( a_i(ii,ij,1) > epsi10 .AND. ht_i(ii,ij,1) < hiclim ) THEN 514 454 a_i(ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / hiclim 515 455 ht_i(ii,ij,1) = hiclim 516 v_i(ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) !clem-useless517 456 ENDIF 518 457 END DO !ji … … 542 481 CALL wrk_dealloc( jpi,jpj, zremap_flag ) ! integer 543 482 CALL wrk_dealloc( jpi,jpj,jpl-1, zdonor ) ! integer 544 CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_ o, dummy_es )483 CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es ) 545 484 CALL wrk_dealloc( jpi,jpj,jpl-1, zdaice, zdvice ) 546 485 CALL wrk_dealloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 ) … … 647 586 REAL(wp) :: zdo_aice ! ice age times volume transferred 648 587 REAL(wp) :: zdaTsf ! aicen*Tsfcn transferred 649 REAL(wp) :: zindsn ! snow or not650 REAL(wp) :: zindb ! ice or not651 588 652 589 INTEGER, POINTER, DIMENSION(:) :: nind_i, nind_j ! compressed indices for i/j directions … … 775 712 776 713 jl1 = zdonor(ii,ij,jl) 777 zindb= MAX( 0.0 , SIGN( 1.0 , v_i(ii,ij,jl1) - epsi10 ) )778 zworka(ii,ij) = zdvice(ii,ij,jl) / MAX(v_i(ii,ij,jl1),epsi10) * zindb714 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) * rswitch 779 716 IF( jl1 == jl) THEN ; jl2 = jl1+1 780 717 ELSE ; jl2 = jl … … 799 736 !-------------- 800 737 801 zdvsnow 738 zdvsnow = v_s(ii,ij,jl1) * zworka(ii,ij) 802 739 v_s(ii,ij,jl1) = v_s(ii,ij,jl1) - zdvsnow 803 740 v_s(ii,ij,jl2) = v_s(ii,ij,jl2) + zdvsnow … … 807 744 !-------------------- 808 745 809 zdesnow 746 zdesnow = e_s(ii,ij,1,jl1) * zworka(ii,ij) 810 747 e_s(ii,ij,1,jl1) = e_s(ii,ij,1,jl1) - zdesnow 811 748 e_s(ii,ij,1,jl2) = e_s(ii,ij,1,jl2) + zdesnow … … 815 752 !-------------- 816 753 817 zdo_aice 754 zdo_aice = oa_i(ii,ij,jl1) * zdaice(ii,ij,jl) 818 755 oa_i(ii,ij,jl1) = oa_i(ii,ij,jl1) - zdo_aice 819 756 oa_i(ii,ij,jl2) = oa_i(ii,ij,jl2) + zdo_aice … … 823 760 !-------------- 824 761 825 zdsm_vice 762 zdsm_vice = smv_i(ii,ij,jl1) * zworka(ii,ij) 826 763 smv_i(ii,ij,jl1) = smv_i(ii,ij,jl1) - zdsm_vice 827 764 smv_i(ii,ij,jl2) = smv_i(ii,ij,jl2) + zdsm_vice … … 831 768 !--------------------- 832 769 833 zdaTsf 770 zdaTsf = t_su(ii,ij,jl1) * zdaice(ii,ij,jl) 834 771 zaTsfn(ii,ij,jl1) = zaTsfn(ii,ij,jl1) - zdaTsf 835 772 zaTsfn(ii,ij,jl2) = zaTsfn(ii,ij,jl2) + zdaTsf … … 872 809 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl) 873 810 t_su(ji,jj,jl) = zaTsfn(ji,jj,jl) / a_i(ji,jj,jl) 874 zindsn= 1.0 - MAX(0.0,SIGN(1.0,-v_s(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes811 rswitch = 1.0 - MAX(0.0,SIGN(1.0,-v_s(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes 875 812 ELSE 876 813 ht_i(ji,jj,jl) = 0._wp … … 888 825 889 826 890 SUBROUTINE lim_itd_th_reb( klbnd, kubnd , ntyp)827 SUBROUTINE lim_itd_th_reb( klbnd, kubnd ) 891 828 !!------------------------------------------------------------------ 892 829 !! *** ROUTINE lim_itd_th_reb *** … … 898 835 INTEGER , INTENT (in) :: klbnd ! Start thickness category index point 899 836 INTEGER , INTENT (in) :: kubnd ! End point on which the the computation is applied 900 INTEGER , INTENT (in) :: ntyp ! number of the ice type involved in the rebinning process901 837 ! 902 838 INTEGER :: ji,jj, jl ! dummy loop indices … … 910 846 REAL(wp), POINTER, DIMENSION(:,:) :: vt_s_init, vt_s_final ! snow volume summed over categories 911 847 !!------------------------------------------------------------------ 848 !! clem 2014/04: be carefull, rebining does not conserve salt(maybe?) => the difference is taken into account in limupdate 912 849 913 850 CALL wrk_alloc( jpi,jpj,jpl, zdonor ) ! interger … … 937 874 938 875 !------------------------------------------------------------------------------ 939 ! 2) Make sure thickness of cat klbnd is at least hi_max _typ(klbnd)876 ! 2) Make sure thickness of cat klbnd is at least hi_max(klbnd) 940 877 !------------------------------------------------------------------------------ 941 878 DO jj = 1, jpj 942 879 DO ji = 1, jpi 943 880 IF( a_i(ji,jj,klbnd) > epsi10 ) THEN 944 IF( ht_i(ji,jj,klbnd) <= hi_max _typ(0,ntyp) .AND. hi_max_typ(0,ntyp) > 0._wp ) THEN945 a_i(ji,jj,klbnd) = v_i(ji,jj,klbnd) / hi_max _typ(0,ntyp)946 ht_i(ji,jj,klbnd) = hi_max _typ(0,ntyp)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) 947 884 ENDIF 948 885 ENDIF … … 1015 952 1016 953 !clem-change 954 DO jj = 1, jpj 955 DO ji = 1, jpi 956 IF( a_i(ji,jj,jl+1) > epsi10 .AND. ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN 957 ! 958 zshiftflag = 1 959 zdonor(ji,jj,jl) = jl + 1 960 zdaice(ji,jj,jl) = a_i(ji,jj,jl+1) 961 zdvice(ji,jj,jl) = v_i(ji,jj,jl+1) 962 ENDIF 963 END DO ! ji 964 END DO ! jj 965 966 IF(lk_mpp) CALL mpp_max( zshiftflag ) 967 968 IF( zshiftflag == 1 ) THEN ! Shift ice between categories 969 CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice ) 970 ! Reset shift parameters 971 zdonor(:,:,jl) = 0 972 zdaice(:,:,jl) = 0._wp 973 zdvice(:,:,jl) = 0._wp 974 ENDIF 975 !clem-change 976 977 ! ! clem-change begin: why not doing that? 1017 978 ! DO jj = 1, jpj 1018 979 ! DO ji = 1, jpi 1019 ! IF( a_i(ji,jj,jl+1) > epsi10 .AND. & 1020 ! ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN 1021 ! ! 1022 ! zshiftflag = 1 1023 ! zdonor(ji,jj,jl) = jl + 1 1024 ! zdaice(ji,jj,jl) = a_i(ji,jj,jl+1) 1025 ! zdvice(ji,jj,jl) = v_i(ji,jj,jl+1) 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) 1026 983 ! ENDIF 1027 984 ! END DO ! ji 1028 985 ! END DO ! jj 1029 !1030 ! IF(lk_mpp) CALL mpp_max( zshiftflag )1031 !1032 ! IF( zshiftflag == 1 ) THEN ! Shift ice between categories1033 ! CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice )1034 ! ! Reset shift parameters1035 ! zdonor(:,:,jl) = 01036 ! zdaice(:,:,jl) = 0._wp1037 ! zdvice(:,:,jl) = 0._wp1038 ! ENDIF1039 !clem-change1040 1041 ! clem-change begin: why not doing that?1042 DO jj = 1, jpj1043 DO ji = 1, jpi1044 IF( a_i(ji,jj,jl+1) > epsi10 .AND. &1045 ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN1046 ht_i(ji,jj,jl+1) = hi_max(jl) + epsi101047 a_i (ji,jj,jl+1) = v_i(ji,jj,jl+1) / ht_i(ji,jj,jl+1)1048 ENDIF1049 END DO ! ji1050 END DO ! jj1051 986 ! clem-change end 1052 987
Note: See TracChangeset
for help on using the changeset viewer.