Changeset 5837 for branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
- Timestamp:
- 2015-10-26T15:59:39+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO14.4_OBS_GENERAL_VINTERP/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
r4333 r5837 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 ! Timing32 USE limcons ! conservation tests 37 33 38 34 IMPLICIT NONE 39 35 PRIVATE 40 36 41 PUBLIC lim_itd_th ! called by ice_stp42 37 PUBLIC lim_itd_th_rem 43 38 PUBLIC lim_itd_th_reb 44 PUBLIC lim_itd_fitline45 PUBLIC lim_itd_shiftice46 47 REAL(wp) :: epsi10 = 1.e-10_wp !48 REAL(wp) :: epsi06 = 1.e-6_wp !49 39 50 40 !!---------------------------------------------------------------------- … … 55 45 CONTAINS 56 46 57 SUBROUTINE lim_itd_th( kt ) 58 !!------------------------------------------------------------------ 59 !! *** ROUTINE lim_itd_th *** 60 !! 61 !! ** Purpose : computes the thermodynamics of ice thickness distribution 62 !! 63 !! ** Method : 64 !!------------------------------------------------------------------ 65 INTEGER, INTENT(in) :: kt ! time step index 66 ! 67 INTEGER :: jl, ja, jm, jbnd1, jbnd2 ! ice types dummy loop index 68 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) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 70 !!------------------------------------------------------------------ 71 IF( nn_timing == 1 ) CALL timing_start('limitd_th') 72 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 ! ------------------------------- 83 84 IF( kt == nit000 .AND. lwp ) THEN 85 WRITE(numout,*) 86 WRITE(numout,*) 'lim_itd_th : Thermodynamics of the ice thickness distribution' 87 WRITE(numout,*) '~~~~~~~~~~~' 88 ENDIF 89 90 !------------------------------------------------------------------------------| 91 ! 1) Transport of ice between thickness categories. | 92 !------------------------------------------------------------------------------| 93 ! Given thermodynamic growth rates, transport ice between 94 ! 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 100 ! 101 CALL lim_var_glo2eqv ! only for info 102 CALL lim_var_agg(1) 103 104 !------------------------------------------------------------------------------| 105 ! 3) Add frazil ice growing in leads. 106 !------------------------------------------------------------------------------| 107 108 CALL lim_thd_lac 109 CALL lim_var_glo2eqv ! only for info 110 111 IF(ln_ctl) THEN ! Control print 112 CALL prt_ctl_info(' ') 113 CALL prt_ctl_info(' - Cell values : ') 114 CALL prt_ctl_info(' ~~~~~~~~~~~~~ ') 115 CALL prt_ctl(tab2d_1=area , clinfo1=' lim_itd_th : cell area :') 116 CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_th : at_i :') 117 CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_th : vt_i :') 118 CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_th : vt_s :') 119 DO jl = 1, jpl 120 CALL prt_ctl_info(' ') 121 CALL prt_ctl_info(' - Category : ', ivar1=jl) 122 CALL prt_ctl_info(' ~~~~~~~~~~') 123 CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_itd_th : a_i : ') 124 CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_itd_th : ht_i : ') 125 CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_itd_th : ht_s : ') 126 CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_itd_th : v_i : ') 127 CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_itd_th : v_s : ') 128 CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_itd_th : e_s : ') 129 CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_itd_th : t_su : ') 130 CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_itd_th : t_snow : ') 131 CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_itd_th : sm_i : ') 132 CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_itd_th : smv_i : ') 133 DO ja = 1, nlay_i 134 CALL prt_ctl_info(' ') 135 CALL prt_ctl_info(' - Layer : ', ivar1=ja) 136 CALL prt_ctl_info(' ~~~~~~~') 137 CALL prt_ctl(tab2d_1=t_i(:,:,ja,jl) , clinfo1= ' lim_itd_th : t_i : ') 138 CALL prt_ctl(tab2d_1=e_i(:,:,ja,jl) , clinfo1= ' lim_itd_th : e_i : ') 139 END DO 140 END DO 141 ENDIF 142 ! 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 ! ------------------------------- 166 ! 167 IF( nn_timing == 1 ) CALL timing_stop('limitd_th') 168 END SUBROUTINE lim_itd_th 169 ! 170 171 SUBROUTINE lim_itd_th_rem( klbnd, kubnd, ntyp, kt ) 47 SUBROUTINE lim_itd_th_rem( klbnd, kubnd, kt ) 172 48 !!------------------------------------------------------------------ 173 49 !! *** ROUTINE lim_itd_th_rem *** … … 182 58 INTEGER , INTENT (in) :: klbnd ! Start thickness category index point 183 59 INTEGER , INTENT (in) :: kubnd ! End point on which the the computation is applied 184 INTEGER , INTENT (in) :: ntyp ! Number of the type used185 60 INTEGER , INTENT (in) :: kt ! Ocean time step 186 61 ! … … 190 65 REAL(wp) :: zx1, zwk1, zdh0, zetamin, zdamax ! local scalars 191 66 REAL(wp) :: zx2, zwk2, zda0, zetamax ! - - 192 REAL(wp) :: zx3 , zareamin, zindb ! - -67 REAL(wp) :: zx3 193 68 CHARACTER (len = 15) :: fieldid 194 69 … … 200 75 REAL(wp), POINTER, DIMENSION(:,:,:) :: hL ! left boundary for the ITD for each thickness 201 76 REAL(wp), POINTER, DIMENSION(:,:,:) :: hR ! left boundary for the ITD for each thickness 202 REAL(wp), POINTER, DIMENSION(:,:,:) :: zht_i_ o! old ice thickness77 REAL(wp), POINTER, DIMENSION(:,:,:) :: zht_i_b ! old ice thickness 203 78 REAL(wp), POINTER, DIMENSION(:,:,:) :: dummy_es 204 79 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdaice, zdvice ! local increment of ice area and volume … … 216 91 !!------------------------------------------------------------------ 217 92 218 CALL wrk_alloc( jpi,jpj, zremap_flag ) ! integer219 CALL wrk_alloc( jpi,jpj,jpl-1, zdonor ) ! integer220 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 ) 221 96 CALL wrk_alloc( jpi,jpj,jpl-1, zdaice, zdvice ) 222 97 CALL wrk_alloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 ) 223 98 CALL wrk_alloc( (jpi+1)*(jpj+1), zvetamin, zvetamax ) 224 CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) ! integer99 CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 225 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 ) 226 227 zareamin = epsi10 !minimum area in thickness categories tolerated by the conceptors of the model228 101 229 102 !!---------------------------------------------------------------------------------------------- … … 247 120 WRITE(numout,*) ' klbnd : ', klbnd 248 121 WRITE(numout,*) ' kubnd : ', kubnd 249 WRITE(numout,*) ' ntyp : ', ntyp250 122 ENDIF 251 123 … … 254 126 DO jj = 1, jpj 255 127 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)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? 261 133 END DO 262 134 END DO … … 277 149 DO jj = 1, jpj 278 150 DO ji = 1, jpi 279 IF ( at_i(ji,jj) .gt. zareamin) THEN151 IF ( at_i(ji,jj) > epsi10 ) THEN 280 152 nbrem = nbrem + 1 281 153 nind_i(nbrem) = ji … … 285 157 zremap_flag(ji,jj) = 0 286 158 ENDIF 287 END DO !ji288 END DO !jj159 END DO 160 END DO 289 161 290 162 !----------------------------------------------------------------------------------------------- … … 292 164 !----------------------------------------------------------------------------------------------- 293 165 !- 4.1 Compute category boundaries 294 ! Tricky trick see limitd_me.F90295 ! will be soon removed, CT296 ! hi_max(kubnd) = 99.297 166 zhbnew(:,:,:) = 0._wp 298 167 … … 302 171 ij = nind_j(ji) 303 172 ! 304 IF ( ( zht_i_o(ii,ij,jl) .GT. epsi10 ) .AND. &305 ( zht_i_o(ii,ij,jl+1) .GT. epsi10 )) THEN173 zhbnew(ii,ij,jl) = hi_max(jl) 174 IF ( a_i_b(ii,ij,jl) > epsi10 .AND. a_i_b(ii,ij,jl+1) > epsi10 ) THEN 306 175 !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 176 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 312 179 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) 313 ELSEIF (zht_i_o(ii,ij,jl+1).gt.epsi10) THEN180 ELSEIF( a_i_b(ii,ij,jl+1) > epsi10) THEN 314 181 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1) 315 ELSE316 zhbnew(ii,ij,jl) = hi_max(jl)317 182 ENDIF 318 183 END DO … … 320 185 !- 4.2 Check that each zhbnew lies between adjacent values of ice thickness 321 186 DO ji = 1, nbrem 322 ! jl, ji323 187 ii = nind_i(ji) 324 188 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 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 329 193 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 194 ELSEIF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) < ( zhbnew(ii,ij,jl) + epsi10 ) ) THEN 333 195 zremap_flag(ii,ij) = 0 334 196 ENDIF 335 197 336 198 !- 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 348 END DO !jl 199 IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0 200 IF( zhbnew(ii,ij,jl) > hi_max(jl+1) ) zremap_flag(ii,ij) = 0 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 349 208 350 209 !----------------------------------------------------------------------------------------------- … … 354 213 DO jj = 1, jpj 355 214 DO ji = 1, jpi 356 IF 215 IF( zremap_flag(ji,jj) == 1 ) THEN 357 216 nbrem = nbrem + 1 358 217 nind_i(nbrem) = ji 359 218 nind_j(nbrem) = jj 360 219 ENDIF 361 END DO !ji362 END DO !jj220 END DO 221 END DO 363 222 364 223 !----------------------------------------------------------------------------------------------- … … 367 226 DO jj = 1, jpj 368 227 DO ji = 1, jpi 369 zhb0(ji,jj) = hi_max_typ(0,ntyp) ! 0eme 370 zhb1(ji,jj) = hi_max_typ(1,ntyp) ! 1er 371 372 zhbnew(ji,jj,klbnd-1) = 0._wp 228 zhb0(ji,jj) = hi_max(0) 229 zhb1(ji,jj) = hi_max(1) 373 230 374 231 IF( a_i(ji,jj,kubnd) > epsi10 ) THEN 375 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) ) 376 233 ELSE 377 zhbnew(ji,jj,kubnd) = hi_max(kubnd) 378 !!? clem bug: since hi_max(jpl)=99, this limit is very high 379 !!? but I think it is erased in fitline subroutine 380 ENDIF 381 382 IF( zhbnew(ji,jj,kubnd) < hi_max(kubnd-1) ) zhbnew(ji,jj,kubnd) = hi_max(kubnd-1) 383 384 END DO !jj 385 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 386 248 387 249 !----------------------------------------------------------------------------------------------- … … 389 251 !----------------------------------------------------------------------------------------------- 390 252 !- 7.1 g(h) for category 1 at start of time step 391 CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_o(:,:,klbnd), & 392 & g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd), & 253 CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd), g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd), & 393 254 & hR(:,:,klbnd), zremap_flag ) 394 255 … … 398 259 ij = nind_j(ji) 399 260 400 !ji401 IF (a_i(ii,ij,klbnd) .gt. epsi10) THEN 261 IF( a_i(ii,ij,klbnd) > epsi10 ) THEN 262 402 263 zdh0 = zdhice(ii,ij,klbnd) !decrease of ice thickness in the lower category 403 ! ji, a_i > epsi10 404 IF (zdh0 .lt. 0.0) THEN !remove area from category 1 405 ! ji, a_i > epsi10; zdh0 < 0 406 zdh0 = MIN(-zdh0,hi_max(klbnd)) 407 264 265 IF( zdh0 < 0.0 ) THEN !remove area from category 1 266 zdh0 = MIN( -zdh0, hi_max(klbnd) ) 408 267 !Integrate g(1) from 0 to dh0 to estimate area melted 409 zetamax = MIN(zdh0,hR(ii,ij,klbnd)) - hL(ii,ij,klbnd) 410 IF (zetamax.gt.0.0) THEN 411 zx1 = zetamax 412 zx2 = 0.5 * zetamax*zetamax 413 zda0 = g1(ii,ij,klbnd) * zx2 + g0(ii,ij,klbnd) * zx1 !ice area removed 414 ! Constrain new thickness <= ht_i 415 zdamax = a_i(ii,ij,klbnd) * & 416 (1.0 - ht_i(ii,ij,klbnd)/zht_i_o(ii,ij,klbnd)) ! zdamax > 0 417 !ice area lost due to melting of thin ice 418 zda0 = MIN(zda0, zdamax) 419 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) 420 277 ! Remove area, conserving volume 421 ht_i(ii,ij,klbnd) = ht_i(ii,ij,klbnd) & 422 * 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 ) 423 279 a_i(ii,ij,klbnd) = a_i(ii,ij,klbnd) - zda0 424 v_i(ii,ij,klbnd) = a_i(ii,ij,klbnd)*ht_i(ii,ij,klbnd) ! clem-useless ? 425 ENDIF ! zetamax > 0 426 ! ji, a_i > epsi10 427 428 ELSE ! if ice accretion 429 ! ji, a_i > epsi10; zdh0 > 0 430 IF ( ntyp .EQ. 1 ) zhbnew(ii,ij,klbnd-1) = MIN(zdh0,hi_max(klbnd)) 431 ! zhbnew was 0, and is shifted to the right to account for thin ice 432 ! growth in openwater (F0 = f1) 433 IF ( ntyp .NE. 1 ) zhbnew(ii,ij,0) = 0 434 ! in other types there is 435 ! no open water growth (F0 = 0) 436 ENDIF ! zdh0 437 438 ! a_i > epsi10 439 ENDIF ! a_i > epsi10 440 441 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 442 291 443 292 !- 7.3 g(h) for each thickness category 444 293 DO jl = klbnd, kubnd 445 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) 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 ) 448 296 END DO 449 297 … … 465 313 ij = nind_j(ji) 466 314 467 IF (zhbnew(ii,ij,jl) .gt. hi_max(jl)) THEN ! transfer from jl to jl+1 468 315 IF (zhbnew(ii,ij,jl) > hi_max(jl)) THEN ! transfer from jl to jl+1 469 316 ! left and right integration limits in eta space 470 zvetamin(ji) = MAX( hi_max(jl), hL(ii,ij,jl)) - hL(ii,ij,jl)471 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) 472 319 zdonor(ii,ij,jl) = jl 473 320 474 ELSE ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl 475 321 ELSE ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl 476 322 ! left and right integration limits in eta space 477 323 zvetamin(ji) = 0.0 478 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) 479 325 zdonor(ii,ij,jl) = jl + 1 480 326 481 ENDIF ! zhbnew(jl) > hi_max(jl)482 483 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 484 330 zetamin = zvetamin(ji) 485 331 486 332 zx1 = zetamax - zetamin 487 zwk1 = zetamin *zetamin488 zwk2 = zetamax *zetamax489 zx2 = 0.5 * ( zwk2 - zwk1)333 zwk1 = zetamin * zetamin 334 zwk2 = zetamax * zetamax 335 zx2 = 0.5 * ( zwk2 - zwk1 ) 490 336 zwk1 = zwk1 * zetamin 491 337 zwk2 = zwk2 * zetamax 492 zx3 = 1.0 /3.0 * (zwk2 - zwk1)338 zx3 = 1.0 / 3.0 * ( zwk2 - zwk1 ) 493 339 nd = zdonor(ii,ij,jl) 494 340 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) 497 498 END DO ! ji 499 END DO ! jl klbnd -> kubnd - 1 341 zdvice(ii,ij,jl) = g1(ii,ij,nd)*zx3 + g0(ii,ij,nd)*zx2 + zdaice(ii,ij,jl)*hL(ii,ij,nd) 342 343 END DO 344 END DO 500 345 501 346 !!---------------------------------------------------------------------------------------------- … … 511 356 ii = nind_i(ji) 512 357 ij = nind_j(ji) 513 IF ( ( a_i(ii,ij,1) > epsi10 ) .AND. ( ht_i(ii,ij,1) < hiclim ) ) THEN 514 a_i(ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / hiclim 515 ht_i(ii,ij,1) = hiclim 516 v_i(ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) !clem-useless 358 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 517 361 ENDIF 518 END DO !ji362 END DO 519 363 520 364 !!---------------------------------------------------------------------------------------------- … … 540 384 ENDIF 541 385 542 CALL wrk_dealloc( jpi,jpj, zremap_flag ) ! integer543 CALL wrk_dealloc( jpi,jpj,jpl-1, zdonor ) ! integer544 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 ) 545 389 CALL wrk_dealloc( jpi,jpj,jpl-1, zdaice, zdvice ) 546 390 CALL wrk_dealloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 ) 547 391 CALL wrk_dealloc( (jpi+1)*(jpj+1), zvetamin, zvetamax ) 548 CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) ! integer392 CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 549 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 ) 550 394 … … 552 396 553 397 554 SUBROUTINE lim_itd_fitline( num_cat, HbL, Hbr, hice, & 555 & g0, g1, hL, hR, zremap_flag ) 398 SUBROUTINE lim_itd_fitline( num_cat, HbL, Hbr, hice, g0, g1, hL, hR, zremap_flag ) 556 399 !!------------------------------------------------------------------ 557 400 !! *** ROUTINE lim_itd_fitline *** … … 572 415 INTEGER , DIMENSION(jpi,jpj), INTENT(in ) :: zremap_flag ! 573 416 ! 574 INTEGER :: ji,jj! horizontal indices417 INTEGER :: ji,jj ! horizontal indices 575 418 REAL(wp) :: zh13 ! HbL + 1/3 * (HbR - HbL) 576 419 REAL(wp) :: zh23 ! HbL + 2/3 * (HbR - HbL) … … 579 422 !!------------------------------------------------------------------ 580 423 ! 581 !582 424 DO jj = 1, jpj 583 425 DO ji = 1, jpi 584 426 ! 585 427 IF( zremap_flag(ji,jj) == 1 .AND. a_i(ji,jj,num_cat) > epsi10 & 586 & .AND. hice(ji,jj) > 0._wp )THEN428 & .AND. hice(ji,jj) > 0._wp ) THEN 587 429 588 430 ! Initialize hL and hR 589 590 431 hL(ji,jj) = HbL(ji,jj) 591 432 hR(ji,jj) = HbR(ji,jj) 592 433 593 434 ! Change hL or hR if hice falls outside central third of range 594 595 zh13 = 1.0/3.0 * (2.0*hL(ji,jj) + hR(ji,jj)) 596 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) ) 597 437 598 438 IF ( hice(ji,jj) < zh13 ) THEN ; hR(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hL(ji,jj) … … 601 441 602 442 ! Compute coefficients of g(eta) = g0 + g1*eta 603 604 443 zdhr = 1._wp / (hR(ji,jj) - hL(ji,jj)) 605 444 zwk1 = 6._wp * a_i(ji,jj,num_cat) * zdhr 606 445 zwk2 = ( hice(ji,jj) - hL(ji,jj) ) * zdhr 607 g0(ji,jj) = zwk1 * ( 2._wp /3._wp - zwk2 )608 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 ) 609 448 ! 610 ELSE 449 ELSE ! remap_flag = .false. or a_i < epsi10 611 450 hL(ji,jj) = 0._wp 612 451 hR(ji,jj) = 0._wp 613 452 g0(ji,jj) = 0._wp 614 453 g1(ji,jj) = 0._wp 615 ENDIF ! a_i > epsi10454 ENDIF 616 455 ! 617 456 END DO … … 637 476 638 477 INTEGER :: ji, jj, jl, jl2, jl1, jk ! dummy loop indices 639 INTEGER :: ii, ij ! indices when changing from 2D-1D is done478 INTEGER :: ii, ij ! indices when changing from 2D-1D is done 640 479 641 480 REAL(wp), POINTER, DIMENSION(:,:,:) :: zaTsfn … … 647 486 REAL(wp) :: zdo_aice ! ice age times volume transferred 648 487 REAL(wp) :: zdaTsf ! aicen*Tsfcn transferred 649 REAL(wp) :: zindsn ! snow or not650 REAL(wp) :: zindb ! ice or not651 488 652 489 INTEGER, POINTER, DIMENSION(:) :: nind_i, nind_j ! compressed indices for i/j directions 653 490 654 INTEGER :: nbrem ! number of cells with ice to transfer 655 656 LOGICAL :: zdaice_negative ! true if daice < -puny 657 LOGICAL :: zdvice_negative ! true if dvice < -puny 658 LOGICAL :: zdaice_greater_aicen ! true if daice > aicen 659 LOGICAL :: zdvice_greater_vicen ! true if dvice > vicen 491 INTEGER :: nbrem ! number of cells with ice to transfer 660 492 !!------------------------------------------------------------------ 661 493 662 494 CALL wrk_alloc( jpi,jpj,jpl, zaTsfn ) 663 495 CALL wrk_alloc( jpi,jpj, zworka ) 664 CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) ! integer496 CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 665 497 666 498 !---------------------------------------------------------------------------------------------- … … 669 501 670 502 DO jl = klbnd, kubnd 671 zaTsfn(:,:,jl) = a_i(:,:,jl)*t_su(:,:,jl) 672 END DO 673 674 !---------------------------------------------------------------------------------------------- 675 ! 2) Check for daice or dvice out of range, allowing for roundoff error 676 !---------------------------------------------------------------------------------------------- 677 ! Note: zdaice < 0 or zdvice < 0 usually happens when category jl 678 ! has a small area, with h(n) very close to a boundary. Then 679 ! the coefficients of g(h) are large, and the computed daice and 680 ! dvice can be in error. If this happens, it is best to transfer 681 ! either the entire category or nothing at all, depending on which 682 ! side of the boundary hice(n) lies. 683 !----------------------------------------------------------------- 684 DO jl = klbnd, kubnd-1 685 686 zdaice_negative = .false. 687 zdvice_negative = .false. 688 zdaice_greater_aicen = .false. 689 zdvice_greater_vicen = .false. 690 691 DO jj = 1, jpj 692 DO ji = 1, jpi 693 694 IF (zdonor(ji,jj,jl) .GT. 0) THEN 695 jl1 = zdonor(ji,jj,jl) 696 697 IF (zdaice(ji,jj,jl) .LT. 0.0) THEN 698 IF (zdaice(ji,jj,jl) .GT. -epsi10) THEN 699 IF ( ( jl1.EQ.jl .AND. ht_i(ji,jj,jl1) .GT. hi_max(jl) ) & 700 .OR. & 701 ( jl1.EQ.jl+1 .AND. ht_i(ji,jj,jl1) .LE. hi_max(jl) ) & 702 ) THEN 703 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category 704 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 705 ELSE 706 zdaice(ji,jj,jl) = 0.0 ! shift no ice 707 zdvice(ji,jj,jl) = 0.0 708 ENDIF 709 ELSE 710 zdaice_negative = .true. 711 ENDIF 712 ENDIF 713 714 IF (zdvice(ji,jj,jl) .LT. 0.0) THEN 715 IF (zdvice(ji,jj,jl) .GT. -epsi10 ) THEN 716 IF ( ( jl1.EQ.jl .AND. ht_i(ji,jj,jl1).GT.hi_max(jl) ) & 717 .OR. & 718 ( jl1.EQ.jl+1 .AND. ht_i(ji,jj,jl1) .LE. hi_max(jl) ) & 719 ) THEN 720 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category 721 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 722 ELSE 723 zdaice(ji,jj,jl) = 0.0 ! shift no ice 724 zdvice(ji,jj,jl) = 0.0 725 ENDIF 726 ELSE 727 zdvice_negative = .true. 728 ENDIF 729 ENDIF 730 731 ! If daice is close to aicen, set daice = aicen. 732 IF (zdaice(ji,jj,jl) .GT. a_i(ji,jj,jl1) - epsi10 ) THEN 733 IF (zdaice(ji,jj,jl) .LT. a_i(ji,jj,jl1)+epsi10) THEN 734 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 735 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 736 ELSE 737 zdaice_greater_aicen = .true. 738 ENDIF 739 ENDIF 740 741 IF (zdvice(ji,jj,jl) .GT. v_i(ji,jj,jl1)-epsi10) THEN 742 IF (zdvice(ji,jj,jl) .LT. v_i(ji,jj,jl1)+epsi10) THEN 743 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 744 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 745 ELSE 746 zdvice_greater_vicen = .true. 747 ENDIF 748 ENDIF 749 750 ENDIF ! donor > 0 751 END DO ! i 752 END DO ! j 753 754 END DO !jl 503 zaTsfn(:,:,jl) = a_i(:,:,jl) * t_su(:,:,jl) 504 END DO 755 505 756 506 !------------------------------------------------------------------------------- 757 ! 3) Transfer volume and energy between categories507 ! 2) Transfer volume and energy between categories 758 508 !------------------------------------------------------------------------------- 759 509 … … 762 512 DO jj = 1, jpj 763 513 DO ji = 1, jpi 764 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 765 515 nbrem = nbrem + 1 766 516 nind_i(nbrem) = ji 767 517 nind_j(nbrem) = jj 768 ENDIF ! tmask518 ENDIF 769 519 END DO 770 520 END DO … … 775 525 776 526 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) * 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 779 529 IF( jl1 == jl) THEN ; jl2 = jl1+1 780 ELSE 530 ELSE ; jl2 = jl 781 531 ENDIF 782 532 … … 784 534 ! Ice areas 785 535 !-------------- 786 787 536 a_i(ii,ij,jl1) = a_i(ii,ij,jl1) - zdaice(ii,ij,jl) 788 537 a_i(ii,ij,jl2) = a_i(ii,ij,jl2) + zdaice(ii,ij,jl) … … 791 540 ! Ice volumes 792 541 !-------------- 793 794 542 v_i(ii,ij,jl1) = v_i(ii,ij,jl1) - zdvice(ii,ij,jl) 795 543 v_i(ii,ij,jl2) = v_i(ii,ij,jl2) + zdvice(ii,ij,jl) … … 798 546 ! Snow volumes 799 547 !-------------- 800 801 zdvsnow = v_s(ii,ij,jl1) * zworka(ii,ij) 548 zdvsnow = v_s(ii,ij,jl1) * zworka(ii,ij) 802 549 v_s(ii,ij,jl1) = v_s(ii,ij,jl1) - zdvsnow 803 550 v_s(ii,ij,jl2) = v_s(ii,ij,jl2) + zdvsnow … … 806 553 ! Snow heat content 807 554 !-------------------- 808 809 zdesnow = e_s(ii,ij,1,jl1) * zworka(ii,ij) 555 zdesnow = e_s(ii,ij,1,jl1) * zworka(ii,ij) 810 556 e_s(ii,ij,1,jl1) = e_s(ii,ij,1,jl1) - zdesnow 811 557 e_s(ii,ij,1,jl2) = e_s(ii,ij,1,jl2) + zdesnow … … 814 560 ! Ice age 815 561 !-------------- 816 817 zdo_aice = oa_i(ii,ij,jl1) * zdaice(ii,ij,jl) 562 zdo_aice = oa_i(ii,ij,jl1) * zdaice(ii,ij,jl) 818 563 oa_i(ii,ij,jl1) = oa_i(ii,ij,jl1) - zdo_aice 819 564 oa_i(ii,ij,jl2) = oa_i(ii,ij,jl2) + zdo_aice … … 822 567 ! Ice salinity 823 568 !-------------- 824 825 zdsm_vice = smv_i(ii,ij,jl1) * zworka(ii,ij) 569 zdsm_vice = smv_i(ii,ij,jl1) * zworka(ii,ij) 826 570 smv_i(ii,ij,jl1) = smv_i(ii,ij,jl1) - zdsm_vice 827 571 smv_i(ii,ij,jl2) = smv_i(ii,ij,jl2) + zdsm_vice … … 830 574 ! Surface temperature 831 575 !--------------------- 832 833 zdaTsf = t_su(ii,ij,jl1) * zdaice(ii,ij,jl) 576 zdaTsf = t_su(ii,ij,jl1) * zdaice(ii,ij,jl) 834 577 zaTsfn(ii,ij,jl1) = zaTsfn(ii,ij,jl1) - zdaTsf 835 578 zaTsfn(ii,ij,jl2) = zaTsfn(ii,ij,jl2) + zdaTsf 836 579 837 END DO ! ji580 END DO 838 581 839 582 !------------------ … … 842 585 843 586 DO jk = 1, nlay_i 844 !CDIR NODEP845 587 DO ji = 1, nbrem 846 588 ii = nind_i(ji) … … 848 590 849 591 jl1 = zdonor(ii,ij,jl) 850 IF (jl1 .EQ.jl) THEN592 IF (jl1 == jl) THEN 851 593 jl2 = jl+1 852 594 ELSE ! n1 = n+1 … … 857 599 e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - zdeice 858 600 e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + zdeice 859 END DO ! ji860 END DO ! jk601 END DO 602 END DO 861 603 862 604 END DO ! boundaries, 1 to ncat-1 … … 872 614 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl) 873 615 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 yes875 616 ELSE 876 617 ht_i(ji,jj,jl) = 0._wp 877 t_su(ji,jj,jl) = rt t618 t_su(ji,jj,jl) = rt0 878 619 ENDIF 879 END DO ! ji880 END DO ! jj881 END DO ! jl620 END DO 621 END DO 622 END DO 882 623 ! 883 624 CALL wrk_dealloc( jpi,jpj,jpl, zaTsfn ) 884 625 CALL wrk_dealloc( jpi,jpj, zworka ) 885 CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) ! integer626 CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 886 627 ! 887 628 END SUBROUTINE lim_itd_shiftice 888 629 889 630 890 SUBROUTINE lim_itd_th_reb( klbnd, kubnd , ntyp)631 SUBROUTINE lim_itd_th_reb( klbnd, kubnd ) 891 632 !!------------------------------------------------------------------ 892 633 !! *** ROUTINE lim_itd_th_reb *** … … 898 639 INTEGER , INTENT (in) :: klbnd ! Start thickness category index point 899 640 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 641 ! 902 642 INTEGER :: ji,jj, jl ! dummy loop indices … … 927 667 DO jj = 1, jpj 928 668 DO ji = 1, jpi 929 IF( a_i(ji,jj,jl) > epsi10 ) THEN 930 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 931 ELSE 932 ht_i(ji,jj,jl) = 0._wp 933 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 934 671 END DO 935 672 END DO … … 937 674 938 675 !------------------------------------------------------------------------------ 939 ! 2) Make sure thickness of cat klbnd is at least hi_max_typ(klbnd) 940 !------------------------------------------------------------------------------ 941 DO jj = 1, jpj 942 DO ji = 1, jpi 943 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 ) THEN 945 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) 947 ENDIF 948 ENDIF 949 END DO 950 END DO 951 952 !------------------------------------------------------------------------------ 953 ! 3) If a category thickness is not in bounds, shift the 676 ! 2) If a category thickness is not in bounds, shift the 954 677 ! entire area, volume, and energy to the neighboring category 955 678 !------------------------------------------------------------------------------ … … 980 703 zdonor(ji,jj,jl) = jl 981 704 ! begin TECLIM change 982 !zdaice(ji,jj,jl) = a_i(ji,jj,jl)983 !zdvice(ji,jj,jl) = v_i(ji,jj,jl)984 705 !zdaice(ji,jj,jl) = a_i(ji,jj,jl) * 0.5_wp 985 706 !zdvice(ji,jj,jl) = v_i(ji,jj,jl)-zdaice(ji,jj,jl)*(hi_max(jl)+hi_max(jl-1)) * 0.5_wp 986 707 ! end TECLIM change 987 708 ! clem: how much of a_i you send in cat sup is somewhat arbitrary 988 zdaice(ji,jj,jl) = a_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - hi_max(jl) + epsi 10 ) / ht_i(ji,jj,jl)989 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 ) 990 711 ENDIF 991 END DO ! ji992 END DO ! jj712 END DO 713 END DO 993 714 IF(lk_mpp) CALL mpp_max( zshiftflag ) 994 715 … … 1001 722 ENDIF 1002 723 ! 1003 END DO ! jl724 END DO 1004 725 1005 726 !---------------------------- … … 1014 735 zshiftflag = 0 1015 736 1016 !clem-change1017 ! DO jj = 1, jpj1018 ! DO ji = 1, jpi1019 ! IF( a_i(ji,jj,jl+1) > epsi10 .AND. &1020 ! ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN1021 ! !1022 ! zshiftflag = 11023 ! zdonor(ji,jj,jl) = jl + 11024 ! zdaice(ji,jj,jl) = a_i(ji,jj,jl+1)1025 ! zdvice(ji,jj,jl) = v_i(ji,jj,jl+1)1026 ! ENDIF1027 ! END DO ! ji1028 ! END DO ! jj1029 !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 737 DO jj = 1, jpj 1043 738 DO ji = 1, jpi 1044 IF( a_i(ji,jj,jl+1) > epsi10 .AND. & 1045 ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN 1046 ht_i(ji,jj,jl+1) = hi_max(jl) + epsi10 1047 a_i (ji,jj,jl+1) = v_i(ji,jj,jl+1) / ht_i(ji,jj,jl+1) 739 IF( a_i(ji,jj,jl+1) > epsi10 .AND. ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN 740 ! 741 zshiftflag = 1 742 zdonor(ji,jj,jl) = jl + 1 743 zdaice(ji,jj,jl) = a_i(ji,jj,jl+1) 744 zdvice(ji,jj,jl) = v_i(ji,jj,jl+1) 1048 745 ENDIF 1049 END DO ! ji 1050 END DO ! jj 1051 ! clem-change end 1052 1053 END DO ! jl 746 END DO 747 END DO 748 749 IF(lk_mpp) CALL mpp_max( zshiftflag ) 750 751 IF( zshiftflag == 1 ) THEN ! Shift ice between categories 752 CALL lim_itd_shiftice( klbnd, kubnd, zdonor, zdaice, zdvice ) 753 ! Reset shift parameters 754 zdonor(:,:,jl) = 0 755 zdaice(:,:,jl) = 0._wp 756 zdvice(:,:,jl) = 0._wp 757 ENDIF 758 759 END DO 1054 760 1055 761 !------------------------------------------------------------------------------ 1056 ! 4) Conservation check762 ! 3) Conservation check 1057 763 !------------------------------------------------------------------------------ 1058 764 … … 1067 773 ENDIF 1068 774 ! 1069 CALL wrk_dealloc( jpi,jpj,jpl, zdonor ) ! interger775 CALL wrk_dealloc( jpi,jpj,jpl, zdonor ) 1070 776 CALL wrk_dealloc( jpi,jpj,jpl, zdaice, zdvice ) 1071 777 CALL wrk_dealloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final ) … … 1078 784 !!---------------------------------------------------------------------- 1079 785 CONTAINS 1080 SUBROUTINE lim_itd_th ! Empty routines1081 END SUBROUTINE lim_itd_th1082 SUBROUTINE lim_itd_th_ini1083 END SUBROUTINE lim_itd_th_ini1084 786 SUBROUTINE lim_itd_th_rem 1085 787 END SUBROUTINE lim_itd_th_rem
Note: See TracChangeset
for help on using the changeset viewer.
