- Timestamp:
- 2015-07-10T13:28:53+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4765_CNRS_agrif/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
r4688 r5581 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 !!---------------------------------------------------------------------- … … 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 48 REAL(wp) :: epsi10 = 1.e-10_wp !49 REAL(wp) :: epsi06 = 1.e-6_wp !50 39 51 40 !!---------------------------------------------------------------------- … … 56 45 CONTAINS 57 46 58 SUBROUTINE lim_itd_th( kt ) 59 !!------------------------------------------------------------------ 60 !! *** ROUTINE lim_itd_th *** 61 !! 62 !! ** Purpose : computes the thermodynamics of ice thickness distribution 63 !! 64 !! ** Method : 65 !!------------------------------------------------------------------ 66 INTEGER, INTENT(in) :: kt ! time step index 67 ! 68 INTEGER :: ji,jj, jk, jl, ja, jm, jbnd1, jbnd2 ! ice types dummy loop index 69 ! 70 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 71 !!------------------------------------------------------------------ 72 IF( nn_timing == 1 ) CALL timing_start('limitd_th') 73 74 ! conservation test 75 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_th', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 76 77 IF( kt == nit000 .AND. lwp ) THEN 78 WRITE(numout,*) 79 WRITE(numout,*) 'lim_itd_th : Thermodynamics of the ice thickness distribution' 80 WRITE(numout,*) '~~~~~~~~~~~' 81 ENDIF 82 83 !------------------------------------------------------------------------------| 84 ! 1) Transport of ice between thickness categories. | 85 !------------------------------------------------------------------------------| 86 ! Given thermodynamic growth rates, transport ice between 87 ! thickness categories. 88 DO jm = 1, jpm 89 jbnd1 = ice_cat_bounds(jm,1) 90 jbnd2 = ice_cat_bounds(jm,2) 91 IF( ice_ncat_types(jm) > 1 ) CALL lim_itd_th_rem( jbnd1, jbnd2, jm, kt ) 92 END DO 93 ! 94 CALL lim_var_glo2eqv ! only for info 95 CALL lim_var_agg(1) 96 97 !------------------------------------------------------------------------------| 98 ! 3) Add frazil ice growing in leads. 99 !------------------------------------------------------------------------------| 100 CALL lim_thd_lac 101 CALL lim_var_glo2eqv ! only for info 102 103 IF(ln_ctl) THEN ! Control print 104 CALL prt_ctl_info(' ') 105 CALL prt_ctl_info(' - Cell values : ') 106 CALL prt_ctl_info(' ~~~~~~~~~~~~~ ') 107 CALL prt_ctl(tab2d_1=area , clinfo1=' lim_itd_th : cell area :') 108 CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_th : at_i :') 109 CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_th : vt_i :') 110 CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_th : vt_s :') 111 DO jl = 1, jpl 112 CALL prt_ctl_info(' ') 113 CALL prt_ctl_info(' - Category : ', ivar1=jl) 114 CALL prt_ctl_info(' ~~~~~~~~~~') 115 CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_itd_th : a_i : ') 116 CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_itd_th : ht_i : ') 117 CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_itd_th : ht_s : ') 118 CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_itd_th : v_i : ') 119 CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_itd_th : v_s : ') 120 CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_itd_th : e_s : ') 121 CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_itd_th : t_su : ') 122 CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_itd_th : t_snow : ') 123 CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_itd_th : sm_i : ') 124 CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_itd_th : smv_i : ') 125 DO ja = 1, nlay_i 126 CALL prt_ctl_info(' ') 127 CALL prt_ctl_info(' - Layer : ', ivar1=ja) 128 CALL prt_ctl_info(' ~~~~~~~') 129 CALL prt_ctl(tab2d_1=t_i(:,:,ja,jl) , clinfo1= ' lim_itd_th : t_i : ') 130 CALL prt_ctl(tab2d_1=e_i(:,:,ja,jl) , clinfo1= ' lim_itd_th : e_i : ') 131 END DO 132 END DO 133 ENDIF 134 ! 135 ! conservation test 136 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_th', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 137 ! 138 IF( nn_timing == 1 ) CALL timing_stop('limitd_th') 139 END SUBROUTINE lim_itd_th 140 ! 141 142 SUBROUTINE lim_itd_th_rem( klbnd, kubnd, ntyp, kt ) 47 SUBROUTINE lim_itd_th_rem( klbnd, kubnd, kt ) 143 48 !!------------------------------------------------------------------ 144 49 !! *** ROUTINE lim_itd_th_rem *** … … 153 58 INTEGER , INTENT (in) :: klbnd ! Start thickness category index point 154 59 INTEGER , INTENT (in) :: kubnd ! End point on which the the computation is applied 155 INTEGER , INTENT (in) :: ntyp ! Number of the type used156 60 INTEGER , INTENT (in) :: kt ! Ocean time step 157 61 ! … … 161 65 REAL(wp) :: zx1, zwk1, zdh0, zetamin, zdamax ! local scalars 162 66 REAL(wp) :: zx2, zwk2, zda0, zetamax ! - - 163 REAL(wp) :: zx3 , zareamin, zindb ! - -67 REAL(wp) :: zx3 164 68 CHARACTER (len = 15) :: fieldid 165 69 … … 171 75 REAL(wp), POINTER, DIMENSION(:,:,:) :: hL ! left boundary for the ITD for each thickness 172 76 REAL(wp), POINTER, DIMENSION(:,:,:) :: hR ! left boundary for the ITD for each thickness 173 REAL(wp), POINTER, DIMENSION(:,:,:) :: zht_i_ o! old ice thickness77 REAL(wp), POINTER, DIMENSION(:,:,:) :: zht_i_b ! old ice thickness 174 78 REAL(wp), POINTER, DIMENSION(:,:,:) :: dummy_es 175 79 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdaice, zdvice ! local increment of ice area and volume … … 187 91 !!------------------------------------------------------------------ 188 92 189 CALL wrk_alloc( jpi,jpj, zremap_flag ) ! integer190 CALL wrk_alloc( jpi,jpj,jpl-1, zdonor ) ! integer191 CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_ o, dummy_es )93 CALL wrk_alloc( jpi,jpj, zremap_flag ) 94 CALL wrk_alloc( jpi,jpj,jpl-1, zdonor ) 95 CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es ) 192 96 CALL wrk_alloc( jpi,jpj,jpl-1, zdaice, zdvice ) 193 97 CALL wrk_alloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 ) 194 98 CALL wrk_alloc( (jpi+1)*(jpj+1), zvetamin, zvetamax ) 195 CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) ! integer99 CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 196 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 ) 197 198 zareamin = epsi10 !minimum area in thickness categories tolerated by the conceptors of the model199 101 200 102 !!---------------------------------------------------------------------------------------------- … … 218 120 WRITE(numout,*) ' klbnd : ', klbnd 219 121 WRITE(numout,*) ' kubnd : ', kubnd 220 WRITE(numout,*) ' ntyp : ', ntyp221 122 ENDIF 222 123 … … 225 126 DO jj = 1, jpj 226 127 DO ji = 1, jpi 227 zindb = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i(ji,jj,jl) +epsi10 ) ) !0 if no ice and 1 if yes228 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zindb229 zindb = 1.0 - MAX( 0.0, SIGN( 1.0, - old_a_i(ji,jj,jl) + epsi10) ) !0 if no ice and 1 if yes230 zht_i_ o(ji,jj,jl) = old_v_i(ji,jj,jl) / MAX( old_a_i(ji,jj,jl), epsi10 ) * zindb231 IF( a_i(ji,jj,jl) > epsi10 ) zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_ o(ji,jj,jl)128 rswitch = MAX( 0.0, SIGN( 1.0, a_i(ji,jj,jl) - epsi10 ) ) !0 if no ice and 1 if yes 129 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * rswitch 130 rswitch = MAX( 0.0, SIGN( 1.0, a_i_b(ji,jj,jl) - epsi10) ) 131 zht_i_b(ji,jj,jl) = v_i_b(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) * rswitch 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? 232 133 END DO 233 134 END DO … … 248 149 DO jj = 1, jpj 249 150 DO ji = 1, jpi 250 IF ( at_i(ji,jj) .gt. zareamin) THEN151 IF ( at_i(ji,jj) > epsi10 ) THEN 251 152 nbrem = nbrem + 1 252 153 nind_i(nbrem) = ji … … 256 157 zremap_flag(ji,jj) = 0 257 158 ENDIF 258 END DO !ji259 END DO !jj159 END DO 160 END DO 260 161 261 162 !----------------------------------------------------------------------------------------------- … … 263 164 !----------------------------------------------------------------------------------------------- 264 165 !- 4.1 Compute category boundaries 265 ! Tricky trick see limitd_me.F90266 ! will be soon removed, CT267 ! hi_max(kubnd) = 99.268 166 zhbnew(:,:,:) = 0._wp 269 167 … … 274 172 ! 275 173 zhbnew(ii,ij,jl) = hi_max(jl) 276 IF ( old_a_i(ii,ij,jl) > epsi10 .AND. old_a_i(ii,ij,jl+1) > epsi10 ) THEN174 IF ( a_i_b(ii,ij,jl) > epsi10 .AND. a_i_b(ii,ij,jl+1) > epsi10 ) THEN 277 175 !interpolate between adjacent category growth rates 278 zslope = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_ o(ii,ij,jl+1) - zht_i_o(ii,ij,jl) )279 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_ o(ii,ij,jl) )280 ELSEIF ( old_a_i(ii,ij,jl) > epsi10) THEN176 zslope = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_b(ii,ij,jl+1) - zht_i_b(ii,ij,jl) ) 177 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_b(ii,ij,jl) ) 178 ELSEIF( a_i_b(ii,ij,jl) > epsi10) THEN 281 179 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) 282 ELSEIF ( old_a_i(ii,ij,jl+1) > epsi10) THEN180 ELSEIF( a_i_b(ii,ij,jl+1) > epsi10) THEN 283 181 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1) 284 182 ENDIF … … 289 187 ii = nind_i(ji) 290 188 ij = nind_j(ji) 291 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 292 193 zremap_flag(ii,ij) = 0 293 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 294 195 zremap_flag(ii,ij) = 0 295 196 ENDIF 296 197 297 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 298 200 IF( zhbnew(ii,ij,jl) > hi_max(jl+1) ) zremap_flag(ii,ij) = 0 299 IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0 300 END DO 301 302 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 303 208 304 209 !----------------------------------------------------------------------------------------------- … … 321 226 DO jj = 1, jpj 322 227 DO ji = 1, jpi 323 zhb0(ji,jj) = hi_max_typ(0,ntyp) ! 0eme 324 zhb1(ji,jj) = hi_max_typ(1,ntyp) ! 1er 325 326 zhbnew(ji,jj,klbnd-1) = 0._wp 228 zhb0(ji,jj) = hi_max(0) 229 zhb1(ji,jj) = hi_max(1) 327 230 328 231 IF( a_i(ji,jj,kubnd) > epsi10 ) THEN 329 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) ) 330 233 ELSE 331 zhbnew(ji,jj,kubnd) = hi_max(kubnd) 332 !!? clem bug: since hi_max(jpl)=99, this limit is very high 333 !!? but I think it is erased in fitline subroutine 334 ENDIF 335 336 IF( zhbnew(ji,jj,kubnd) < hi_max(kubnd-1) ) zhbnew(ji,jj,kubnd) = hi_max(kubnd-1) 337 338 END DO !jj 339 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 340 248 341 249 !----------------------------------------------------------------------------------------------- … … 343 251 !----------------------------------------------------------------------------------------------- 344 252 !- 7.1 g(h) for category 1 at start of time step 345 CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_o(:,:,klbnd), & 346 & g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd), & 253 CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd), g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd), & 347 254 & hR(:,:,klbnd), zremap_flag ) 348 255 … … 352 259 ij = nind_j(ji) 353 260 354 !ji355 IF (a_i(ii,ij,klbnd) .gt. epsi10) THEN 261 IF( a_i(ii,ij,klbnd) > epsi10 ) THEN 262 356 263 zdh0 = zdhice(ii,ij,klbnd) !decrease of ice thickness in the lower category 357 ! ji, a_i > epsi10 358 IF (zdh0 .lt. 0.0) THEN !remove area from category 1 359 ! ji, a_i > epsi10; zdh0 < 0 360 zdh0 = MIN(-zdh0,hi_max(klbnd)) 361 264 265 IF( zdh0 < 0.0 ) THEN !remove area from category 1 266 zdh0 = MIN( -zdh0, hi_max(klbnd) ) 362 267 !Integrate g(1) from 0 to dh0 to estimate area melted 363 zetamax = MIN(zdh0,hR(ii,ij,klbnd)) - hL(ii,ij,klbnd) 364 IF (zetamax.gt.0.0) THEN 365 zx1 = zetamax 366 zx2 = 0.5 * zetamax*zetamax 367 zda0 = g1(ii,ij,klbnd) * zx2 + g0(ii,ij,klbnd) * zx1 !ice area removed 368 ! Constrain new thickness <= ht_i 369 zdamax = a_i(ii,ij,klbnd) * & 370 (1.0 - ht_i(ii,ij,klbnd)/zht_i_o(ii,ij,klbnd)) ! zdamax > 0 371 !ice area lost due to melting of thin ice 372 zda0 = MIN(zda0, zdamax) 373 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) 374 277 ! Remove area, conserving volume 375 ht_i(ii,ij,klbnd) = ht_i(ii,ij,klbnd) & 376 * 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 ) 377 279 a_i(ii,ij,klbnd) = a_i(ii,ij,klbnd) - zda0 378 v_i(ii,ij,klbnd) = a_i(ii,ij,klbnd)*ht_i(ii,ij,klbnd) ! clem-useless ? 379 ENDIF ! zetamax > 0 380 ! ji, a_i > epsi10 381 382 ELSE ! if ice accretion 383 ! ji, a_i > epsi10; zdh0 > 0 384 IF ( ntyp .EQ. 1 ) zhbnew(ii,ij,klbnd-1) = MIN(zdh0,hi_max(klbnd)) 385 ! zhbnew was 0, and is shifted to the right to account for thin ice 386 ! growth in openwater (F0 = f1) 387 IF ( ntyp .NE. 1 ) zhbnew(ii,ij,0) = 0 388 ! in other types there is 389 ! no open water growth (F0 = 0) 390 ENDIF ! zdh0 391 392 ! a_i > epsi10 393 ENDIF ! a_i > epsi10 394 395 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 396 291 397 292 !- 7.3 g(h) for each thickness category 398 293 DO jl = klbnd, kubnd 399 CALL lim_itd_fitline( jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl),&400 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 ) 401 296 END DO 402 297 … … 418 313 ij = nind_j(ji) 419 314 420 IF (zhbnew(ii,ij,jl) .gt. hi_max(jl)) THEN ! transfer from jl to jl+1 421 315 IF (zhbnew(ii,ij,jl) > hi_max(jl)) THEN ! transfer from jl to jl+1 422 316 ! left and right integration limits in eta space 423 zvetamin(ji) = MAX( hi_max(jl), hL(ii,ij,jl)) - hL(ii,ij,jl)424 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) 425 319 zdonor(ii,ij,jl) = jl 426 320 427 ELSE ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl 428 321 ELSE ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl 429 322 ! left and right integration limits in eta space 430 323 zvetamin(ji) = 0.0 431 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) 432 325 zdonor(ii,ij,jl) = jl + 1 433 326 434 ENDIF ! zhbnew(jl) > hi_max(jl)435 436 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 437 330 zetamin = zvetamin(ji) 438 331 439 332 zx1 = zetamax - zetamin 440 zwk1 = zetamin *zetamin441 zwk2 = zetamax *zetamax442 zx2 = 0.5 * ( zwk2 - zwk1)333 zwk1 = zetamin * zetamin 334 zwk2 = zetamax * zetamax 335 zx2 = 0.5 * ( zwk2 - zwk1 ) 443 336 zwk1 = zwk1 * zetamin 444 337 zwk2 = zwk2 * zetamax 445 zx3 = 1.0 /3.0 * (zwk2 - zwk1)338 zx3 = 1.0 / 3.0 * ( zwk2 - zwk1 ) 446 339 nd = zdonor(ii,ij,jl) 447 340 zdaice(ii,ij,jl) = g1(ii,ij,nd)*zx2 + g0(ii,ij,nd)*zx1 448 341 zdvice(ii,ij,jl) = g1(ii,ij,nd)*zx3 + g0(ii,ij,nd)*zx2 + zdaice(ii,ij,jl)*hL(ii,ij,nd) 449 342 450 END DO ! ji451 END DO ! jl klbnd -> kubnd - 1343 END DO 344 END DO 452 345 453 346 !!---------------------------------------------------------------------------------------------- … … 463 356 ii = nind_i(ji) 464 357 ij = nind_j(ji) 465 IF ( a_i(ii,ij,1) > epsi10 .AND. ht_i(ii,ij,1) < hiclim) THEN466 a_i (ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / hiclim467 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 468 361 ENDIF 469 END DO !ji362 END DO 470 363 471 364 !!---------------------------------------------------------------------------------------------- … … 491 384 ENDIF 492 385 493 CALL wrk_dealloc( jpi,jpj, zremap_flag ) ! integer494 CALL wrk_dealloc( jpi,jpj,jpl-1, zdonor ) ! integer495 CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_ o, dummy_es )386 CALL wrk_dealloc( jpi,jpj, zremap_flag ) 387 CALL wrk_dealloc( jpi,jpj,jpl-1, zdonor ) 388 CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es ) 496 389 CALL wrk_dealloc( jpi,jpj,jpl-1, zdaice, zdvice ) 497 390 CALL wrk_dealloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 ) 498 391 CALL wrk_dealloc( (jpi+1)*(jpj+1), zvetamin, zvetamax ) 499 CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) ! integer392 CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 500 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 ) 501 394 … … 503 396 504 397 505 SUBROUTINE lim_itd_fitline( num_cat, HbL, Hbr, hice, & 506 & g0, g1, hL, hR, zremap_flag ) 398 SUBROUTINE lim_itd_fitline( num_cat, HbL, Hbr, hice, g0, g1, hL, hR, zremap_flag ) 507 399 !!------------------------------------------------------------------ 508 400 !! *** ROUTINE lim_itd_fitline *** … … 523 415 INTEGER , DIMENSION(jpi,jpj), INTENT(in ) :: zremap_flag ! 524 416 ! 525 INTEGER :: ji,jj! horizontal indices417 INTEGER :: ji,jj ! horizontal indices 526 418 REAL(wp) :: zh13 ! HbL + 1/3 * (HbR - HbL) 527 419 REAL(wp) :: zh23 ! HbL + 2/3 * (HbR - HbL) … … 530 422 !!------------------------------------------------------------------ 531 423 ! 532 !533 424 DO jj = 1, jpj 534 425 DO ji = 1, jpi 535 426 ! 536 427 IF( zremap_flag(ji,jj) == 1 .AND. a_i(ji,jj,num_cat) > epsi10 & 537 & .AND. hice(ji,jj) > 0._wp )THEN428 & .AND. hice(ji,jj) > 0._wp ) THEN 538 429 539 430 ! Initialize hL and hR 540 541 431 hL(ji,jj) = HbL(ji,jj) 542 432 hR(ji,jj) = HbR(ji,jj) 543 433 544 434 ! Change hL or hR if hice falls outside central third of range 545 546 zh13 = 1.0/3.0 * (2.0*hL(ji,jj) + hR(ji,jj)) 547 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) ) 548 437 549 438 IF ( hice(ji,jj) < zh13 ) THEN ; hR(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hL(ji,jj) … … 552 441 553 442 ! Compute coefficients of g(eta) = g0 + g1*eta 554 555 443 zdhr = 1._wp / (hR(ji,jj) - hL(ji,jj)) 556 444 zwk1 = 6._wp * a_i(ji,jj,num_cat) * zdhr 557 445 zwk2 = ( hice(ji,jj) - hL(ji,jj) ) * zdhr 558 g0(ji,jj) = zwk1 * ( 2._wp /3._wp - zwk2 )559 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 ) 560 448 ! 561 ELSE 449 ELSE ! remap_flag = .false. or a_i < epsi10 562 450 hL(ji,jj) = 0._wp 563 451 hR(ji,jj) = 0._wp 564 452 g0(ji,jj) = 0._wp 565 453 g1(ji,jj) = 0._wp 566 ENDIF ! a_i > epsi10454 ENDIF 567 455 ! 568 456 END DO … … 588 476 589 477 INTEGER :: ji, jj, jl, jl2, jl1, jk ! dummy loop indices 590 INTEGER :: ii, ij ! indices when changing from 2D-1D is done478 INTEGER :: ii, ij ! indices when changing from 2D-1D is done 591 479 592 480 REAL(wp), POINTER, DIMENSION(:,:,:) :: zaTsfn … … 598 486 REAL(wp) :: zdo_aice ! ice age times volume transferred 599 487 REAL(wp) :: zdaTsf ! aicen*Tsfcn transferred 600 REAL(wp) :: zindsn ! snow or not601 REAL(wp) :: zindb ! ice or not602 488 603 489 INTEGER, POINTER, DIMENSION(:) :: nind_i, nind_j ! compressed indices for i/j directions 604 490 605 INTEGER :: nbrem ! number of cells with ice to transfer 606 607 LOGICAL :: zdaice_negative ! true if daice < -puny 608 LOGICAL :: zdvice_negative ! true if dvice < -puny 609 LOGICAL :: zdaice_greater_aicen ! true if daice > aicen 610 LOGICAL :: zdvice_greater_vicen ! true if dvice > vicen 491 INTEGER :: nbrem ! number of cells with ice to transfer 611 492 !!------------------------------------------------------------------ 612 493 613 494 CALL wrk_alloc( jpi,jpj,jpl, zaTsfn ) 614 495 CALL wrk_alloc( jpi,jpj, zworka ) 615 CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) ! integer496 CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 616 497 617 498 !---------------------------------------------------------------------------------------------- … … 620 501 621 502 DO jl = klbnd, kubnd 622 zaTsfn(:,:,jl) = a_i(:,:,jl)*t_su(:,:,jl) 623 END DO 624 625 !---------------------------------------------------------------------------------------------- 626 ! 2) Check for daice or dvice out of range, allowing for roundoff error 627 !---------------------------------------------------------------------------------------------- 628 ! Note: zdaice < 0 or zdvice < 0 usually happens when category jl 629 ! has a small area, with h(n) very close to a boundary. Then 630 ! the coefficients of g(h) are large, and the computed daice and 631 ! dvice can be in error. If this happens, it is best to transfer 632 ! either the entire category or nothing at all, depending on which 633 ! side of the boundary hice(n) lies. 634 !----------------------------------------------------------------- 635 DO jl = klbnd, kubnd-1 636 637 zdaice_negative = .false. 638 zdvice_negative = .false. 639 zdaice_greater_aicen = .false. 640 zdvice_greater_vicen = .false. 641 642 DO jj = 1, jpj 643 DO ji = 1, jpi 644 645 IF (zdonor(ji,jj,jl) .GT. 0) THEN 646 jl1 = zdonor(ji,jj,jl) 647 648 IF (zdaice(ji,jj,jl) .LT. 0.0) THEN 649 IF (zdaice(ji,jj,jl) .GT. -epsi10) THEN 650 IF ( ( jl1.EQ.jl .AND. ht_i(ji,jj,jl1) .GT. hi_max(jl) ) & 651 .OR. & 652 ( jl1.EQ.jl+1 .AND. ht_i(ji,jj,jl1) .LE. hi_max(jl) ) & 653 ) THEN 654 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category 655 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 656 ELSE 657 zdaice(ji,jj,jl) = 0.0 ! shift no ice 658 zdvice(ji,jj,jl) = 0.0 659 ENDIF 660 ELSE 661 zdaice_negative = .true. 662 ENDIF 663 ENDIF 664 665 IF (zdvice(ji,jj,jl) .LT. 0.0) THEN 666 IF (zdvice(ji,jj,jl) .GT. -epsi10 ) THEN 667 IF ( ( jl1.EQ.jl .AND. ht_i(ji,jj,jl1).GT.hi_max(jl) ) & 668 .OR. & 669 ( jl1.EQ.jl+1 .AND. ht_i(ji,jj,jl1) .LE. hi_max(jl) ) & 670 ) THEN 671 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category 672 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 673 ELSE 674 zdaice(ji,jj,jl) = 0.0 ! shift no ice 675 zdvice(ji,jj,jl) = 0.0 676 ENDIF 677 ELSE 678 zdvice_negative = .true. 679 ENDIF 680 ENDIF 681 682 ! If daice is close to aicen, set daice = aicen. 683 IF (zdaice(ji,jj,jl) .GT. a_i(ji,jj,jl1) - epsi10 ) THEN 684 IF (zdaice(ji,jj,jl) .LT. a_i(ji,jj,jl1)+epsi10) THEN 685 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 686 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 687 ELSE 688 zdaice_greater_aicen = .true. 689 ENDIF 690 ENDIF 691 692 IF (zdvice(ji,jj,jl) .GT. v_i(ji,jj,jl1)-epsi10) THEN 693 IF (zdvice(ji,jj,jl) .LT. v_i(ji,jj,jl1)+epsi10) THEN 694 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 695 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 696 ELSE 697 zdvice_greater_vicen = .true. 698 ENDIF 699 ENDIF 700 701 ENDIF ! donor > 0 702 END DO ! i 703 END DO ! j 704 705 END DO !jl 503 zaTsfn(:,:,jl) = a_i(:,:,jl) * t_su(:,:,jl) 504 END DO 706 505 707 506 !------------------------------------------------------------------------------- 708 ! 3) Transfer volume and energy between categories507 ! 2) Transfer volume and energy between categories 709 508 !------------------------------------------------------------------------------- 710 509 … … 713 512 DO jj = 1, jpj 714 513 DO ji = 1, jpi 715 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 716 515 nbrem = nbrem + 1 717 516 nind_i(nbrem) = ji 718 517 nind_j(nbrem) = jj 719 ENDIF ! tmask518 ENDIF 720 519 END DO 721 520 END DO … … 726 525 727 526 jl1 = zdonor(ii,ij,jl) 728 zindb = MAX( 0.0 , SIGN( 1.0, v_i(ii,ij,jl1) - epsi10 ) )729 zworka(ii,ij) = zdvice(ii,ij,jl) / MAX(v_i(ii,ij,jl1),epsi10) * zindb527 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 730 529 IF( jl1 == jl) THEN ; jl2 = jl1+1 731 ELSE 530 ELSE ; jl2 = jl 732 531 ENDIF 733 532 … … 735 534 ! Ice areas 736 535 !-------------- 737 738 536 a_i(ii,ij,jl1) = a_i(ii,ij,jl1) - zdaice(ii,ij,jl) 739 537 a_i(ii,ij,jl2) = a_i(ii,ij,jl2) + zdaice(ii,ij,jl) … … 742 540 ! Ice volumes 743 541 !-------------- 744 745 542 v_i(ii,ij,jl1) = v_i(ii,ij,jl1) - zdvice(ii,ij,jl) 746 543 v_i(ii,ij,jl2) = v_i(ii,ij,jl2) + zdvice(ii,ij,jl) … … 749 546 ! Snow volumes 750 547 !-------------- 751 752 548 zdvsnow = v_s(ii,ij,jl1) * zworka(ii,ij) 753 549 v_s(ii,ij,jl1) = v_s(ii,ij,jl1) - zdvsnow … … 757 553 ! Snow heat content 758 554 !-------------------- 759 760 555 zdesnow = e_s(ii,ij,1,jl1) * zworka(ii,ij) 761 556 e_s(ii,ij,1,jl1) = e_s(ii,ij,1,jl1) - zdesnow … … 765 560 ! Ice age 766 561 !-------------- 767 768 562 zdo_aice = oa_i(ii,ij,jl1) * zdaice(ii,ij,jl) 769 563 oa_i(ii,ij,jl1) = oa_i(ii,ij,jl1) - zdo_aice … … 773 567 ! Ice salinity 774 568 !-------------- 775 776 569 zdsm_vice = smv_i(ii,ij,jl1) * zworka(ii,ij) 777 570 smv_i(ii,ij,jl1) = smv_i(ii,ij,jl1) - zdsm_vice … … 781 574 ! Surface temperature 782 575 !--------------------- 783 784 576 zdaTsf = t_su(ii,ij,jl1) * zdaice(ii,ij,jl) 785 577 zaTsfn(ii,ij,jl1) = zaTsfn(ii,ij,jl1) - zdaTsf 786 578 zaTsfn(ii,ij,jl2) = zaTsfn(ii,ij,jl2) + zdaTsf 787 579 788 END DO ! ji580 END DO 789 581 790 582 !------------------ … … 793 585 794 586 DO jk = 1, nlay_i 795 !CDIR NODEP796 587 DO ji = 1, nbrem 797 588 ii = nind_i(ji) … … 799 590 800 591 jl1 = zdonor(ii,ij,jl) 801 IF (jl1 .EQ.jl) THEN592 IF (jl1 == jl) THEN 802 593 jl2 = jl+1 803 594 ELSE ! n1 = n+1 … … 808 599 e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - zdeice 809 600 e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + zdeice 810 END DO ! ji811 END DO ! jk601 END DO 602 END DO 812 603 813 604 END DO ! boundaries, 1 to ncat-1 … … 823 614 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl) 824 615 t_su(ji,jj,jl) = zaTsfn(ji,jj,jl) / a_i(ji,jj,jl) 825 zindsn = 1.0 - MAX(0.0,SIGN(1.0,-v_s(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes826 616 ELSE 827 617 ht_i(ji,jj,jl) = 0._wp 828 t_su(ji,jj,jl) = rt t618 t_su(ji,jj,jl) = rt0 829 619 ENDIF 830 END DO ! ji831 END DO ! jj832 END DO ! jl620 END DO 621 END DO 622 END DO 833 623 ! 834 624 CALL wrk_dealloc( jpi,jpj,jpl, zaTsfn ) 835 625 CALL wrk_dealloc( jpi,jpj, zworka ) 836 CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) ! integer626 CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 837 627 ! 838 628 END SUBROUTINE lim_itd_shiftice 839 629 840 630 841 SUBROUTINE lim_itd_th_reb( klbnd, kubnd , ntyp)631 SUBROUTINE lim_itd_th_reb( klbnd, kubnd ) 842 632 !!------------------------------------------------------------------ 843 633 !! *** ROUTINE lim_itd_th_reb *** … … 849 639 INTEGER , INTENT (in) :: klbnd ! Start thickness category index point 850 640 INTEGER , INTENT (in) :: kubnd ! End point on which the the computation is applied 851 INTEGER , INTENT (in) :: ntyp ! number of the ice type involved in the rebinning process852 641 ! 853 642 INTEGER :: ji,jj, jl ! dummy loop indices … … 861 650 REAL(wp), POINTER, DIMENSION(:,:) :: vt_s_init, vt_s_final ! snow volume summed over categories 862 651 !!------------------------------------------------------------------ 863 !! clem 2014/04: be carefull, rebining does not conserve salt(maybe?) => the difference is taken into account in limupdate864 652 865 653 CALL wrk_alloc( jpi,jpj,jpl, zdonor ) ! interger … … 879 667 DO jj = 1, jpj 880 668 DO ji = 1, jpi 881 IF( a_i(ji,jj,jl) > epsi10 ) THEN 882 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 883 ELSE 884 ht_i(ji,jj,jl) = 0._wp 885 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 886 671 END DO 887 672 END DO … … 889 674 890 675 !------------------------------------------------------------------------------ 891 ! 2) Make sure thickness of cat klbnd is at least hi_max_typ(klbnd) 892 !------------------------------------------------------------------------------ 893 DO jj = 1, jpj 894 DO ji = 1, jpi 895 IF( a_i(ji,jj,klbnd) > epsi10 ) THEN 896 IF( ht_i(ji,jj,klbnd) <= hi_max_typ(0,ntyp) .AND. hi_max_typ(0,ntyp) > 0._wp ) THEN 897 a_i(ji,jj,klbnd) = v_i(ji,jj,klbnd) / hi_max_typ(0,ntyp) 898 ht_i(ji,jj,klbnd) = hi_max_typ(0,ntyp) 899 ENDIF 900 ENDIF 901 END DO 902 END DO 903 904 !------------------------------------------------------------------------------ 905 ! 3) If a category thickness is not in bounds, shift the 676 ! 2) If a category thickness is not in bounds, shift the 906 677 ! entire area, volume, and energy to the neighboring category 907 678 !------------------------------------------------------------------------------ … … 932 703 zdonor(ji,jj,jl) = jl 933 704 ! begin TECLIM change 934 !zdaice(ji,jj,jl) = a_i(ji,jj,jl)935 !zdvice(ji,jj,jl) = v_i(ji,jj,jl)936 705 !zdaice(ji,jj,jl) = a_i(ji,jj,jl) * 0.5_wp 937 706 !zdvice(ji,jj,jl) = v_i(ji,jj,jl)-zdaice(ji,jj,jl)*(hi_max(jl)+hi_max(jl-1)) * 0.5_wp 938 707 ! end TECLIM change 939 708 ! clem: how much of a_i you send in cat sup is somewhat arbitrary 940 zdaice(ji,jj,jl) = a_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - hi_max(jl) + epsi 10 ) / ht_i(ji,jj,jl)941 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 ) 942 711 ENDIF 943 END DO ! ji944 END DO ! jj712 END DO 713 END DO 945 714 IF(lk_mpp) CALL mpp_max( zshiftflag ) 946 715 … … 953 722 ENDIF 954 723 ! 955 END DO ! jl724 END DO 956 725 957 726 !---------------------------- … … 966 735 zshiftflag = 0 967 736 968 !clem-change969 737 DO jj = 1, jpj 970 738 DO ji = 1, jpi … … 976 744 zdvice(ji,jj,jl) = v_i(ji,jj,jl+1) 977 745 ENDIF 978 END DO ! ji979 END DO ! jj746 END DO 747 END DO 980 748 981 749 IF(lk_mpp) CALL mpp_max( zshiftflag ) … … 988 756 zdvice(:,:,jl) = 0._wp 989 757 ENDIF 990 !clem-change 991 992 ! ! clem-change begin: why not doing that? 993 ! DO jj = 1, jpj 994 ! DO ji = 1, jpi 995 ! IF( a_i(ji,jj,jl+1) > epsi10 .AND. ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN 996 ! ht_i(ji,jj,jl+1) = hi_max(jl) + epsi10 997 ! a_i (ji,jj,jl+1) = v_i(ji,jj,jl+1) / ht_i(ji,jj,jl+1) 998 ! ENDIF 999 ! END DO ! ji 1000 ! END DO ! jj 1001 ! clem-change end 1002 1003 END DO ! jl 758 759 END DO 1004 760 1005 761 !------------------------------------------------------------------------------ 1006 ! 4) Conservation check762 ! 3) Conservation check 1007 763 !------------------------------------------------------------------------------ 1008 764 … … 1017 773 ENDIF 1018 774 ! 1019 CALL wrk_dealloc( jpi,jpj,jpl, zdonor ) ! interger775 CALL wrk_dealloc( jpi,jpj,jpl, zdonor ) 1020 776 CALL wrk_dealloc( jpi,jpj,jpl, zdaice, zdvice ) 1021 777 CALL wrk_dealloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final ) … … 1028 784 !!---------------------------------------------------------------------- 1029 785 CONTAINS 1030 SUBROUTINE lim_itd_th ! Empty routines1031 END SUBROUTINE lim_itd_th1032 SUBROUTINE lim_itd_th_ini1033 END SUBROUTINE lim_itd_th_ini1034 786 SUBROUTINE lim_itd_th_rem 1035 787 END SUBROUTINE lim_itd_th_rem
Note: See TracChangeset
for help on using the changeset viewer.
