Changeset 921 for trunk/NEMO/LIM_SRC_3/limitd_th.F90
- Timestamp:
- 2008-05-13T10:28:52+02:00 (16 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/LIM_SRC_3/limitd_th.F90 (modified) (32 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/LIM_SRC_3/limitd_th.F90
r869 r921 28 28 USE prtctl ! Print control 29 29 USE lib_mpp 30 30 31 31 IMPLICIT NONE 32 32 PRIVATE … … 51 51 !!---------------------------------------------------------------------- 52 52 53 !!----------------------------------------------------------------------------------------------54 !!----------------------------------------------------------------------------------------------53 !!---------------------------------------------------------------------------------------------- 54 !!---------------------------------------------------------------------------------------------- 55 55 56 56 CONTAINS 57 57 58 SUBROUTINE lim_itd_th 59 !!------------------------------------------------------------------ 60 !! *** ROUTINE lim_itd_th *** 61 !! ** Purpose : 62 !! This routine computes the thermodynamics of ice thickness 63 !! distribution 64 !! ** Method : 65 !! 66 !! ** Arguments : 67 !! kideb , kiut : Starting and ending points on which the 68 !! the computation is applied 69 !! 70 !! ** Inputs / Ouputs : (global commons) 71 !! 72 !! ** External : 73 !! 74 !! ** References : 75 !! 76 !! ** History : 77 !! (12-2005) Martin Vancoppenolle 78 !! 79 !!------------------------------------------------------------------ 80 !! * Arguments 81 82 !! * Local variables 83 INTEGER :: jl, ja, & ! ice category, layers 84 jm, & ! ice types dummy loop index 85 jbnd1, & 86 jbnd2 87 88 REAL(wp) :: & ! constant values 89 zeps = 1.0e-10, & 90 epsi10 = 1.0e-10 91 92 !!-- End of declarations 93 !!---------------------------------------------------------------------------------------------- 94 95 IF (lwp) THEN 96 WRITE(numout,*) 97 WRITE(numout,*) 'lim_itd_th : Thermodynamics of the ice thickness distribution' 98 WRITE(numout,*) '~~~~~~~~~~~' 99 ENDIF 100 101 !------------------------------------------------------------------------------| 102 ! 1) Transport of ice between thickness categories. | 103 !------------------------------------------------------------------------------| 58 SUBROUTINE lim_itd_th( kt ) 59 !!------------------------------------------------------------------ 60 !! *** ROUTINE lim_itd_th *** 61 !! ** Purpose : 62 !! This routine computes the thermodynamics of ice thickness 63 !! distribution 64 !! ** Method : 65 !! 66 !! ** Arguments : 67 !! kideb , kiut : Starting and ending points on which the 68 !! the computation is applied 69 !! 70 !! ** Inputs / Ouputs : (global commons) 71 !! 72 !! ** External : 73 !! 74 !! ** References : 75 !! 76 !! ** History : 77 !! (12-2005) Martin Vancoppenolle 78 !! 79 !!------------------------------------------------------------------ 80 !! * Arguments 81 INTEGER, INTENT(in) :: kt 82 !! * Local variables 83 INTEGER :: jl, ja, & ! ice category, layers 84 jm, & ! ice types dummy loop index 85 jbnd1, & 86 jbnd2 87 88 REAL(wp) :: & ! constant values 89 zeps = 1.0e-10, & 90 epsi10 = 1.0e-10 91 92 IF( kt == nit000 .AND. lwp ) THEN 93 WRITE(numout,*) 94 WRITE(numout,*) 'lim_itd_th : Thermodynamics of the ice thickness distribution' 95 WRITE(numout,*) '~~~~~~~~~~~' 96 ENDIF 97 98 !------------------------------------------------------------------------------| 99 ! 1) Transport of ice between thickness categories. | 100 !------------------------------------------------------------------------------| 104 101 ! Given thermodynamic growth rates, transport ice between 105 102 ! thickness categories. … … 107 104 jbnd1 = ice_cat_bounds(jm,1) 108 105 jbnd2 = ice_cat_bounds(jm,2) 109 IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_rem( jbnd1, jbnd2, jm)106 IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_rem( jbnd1, jbnd2, jm, kt ) 110 107 END DO 111 108 … … 113 110 CALL lim_var_agg(1) 114 111 115 !------------------------------------------------------------------------------|116 ! 3) Add frazil ice growing in leads.117 !------------------------------------------------------------------------------|112 !------------------------------------------------------------------------------| 113 ! 3) Add frazil ice growing in leads. 114 !------------------------------------------------------------------------------| 118 115 119 116 CALL lim_thd_lac 120 117 CALL lim_var_glo2eqv ! only for info 121 118 122 !----------------------------------------------------------------------------------------123 ! 4) Computation of trend terms and get back to old values124 !----------------------------------------------------------------------------------------119 !---------------------------------------------------------------------------------------- 120 ! 4) Computation of trend terms and get back to old values 121 !---------------------------------------------------------------------------------------- 125 122 126 123 !- Trend terms … … 133 130 d_smv_i_thd(:,:,:) = 0.0 134 131 IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) & 135 d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:)132 d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 136 133 137 134 IF(ln_ctl) THEN ! Control print … … 166 163 END DO 167 164 ENDIF 168 165 169 166 !- Recover Old values 170 167 a_i(:,:,:) = old_a_i (:,:,:) … … 175 172 176 173 IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) & 177 smv_i(:,:,:) = old_smv_i (:,:,:) 178 179 180 END SUBROUTINE lim_itd_th 181 182 !!---------------------------------------------------------------------------------------------- 183 !!---------------------------------------------------------------------------------------------- 184 185 SUBROUTINE lim_itd_th_rem(klbnd,kubnd,ntyp) 186 !!------------------------------------------------------------------ 187 !! *** ROUTINE lim_itd_th_rem *** 188 !! ** Purpose : 189 !! This routine computes the redistribution of ice thickness 190 !! after thermodynamic growth of ice thickness 191 !! 192 !! ** Method : Linear remapping 193 !! 194 !! ** Arguments : 195 !! klbnd, kubnd : Starting and ending category index on which the 196 !! the computation is applied 197 !! 198 !! ** Inputs / Ouputs : (global commons) 199 !! 200 !! ** External : 201 !! 202 !! ** References : W.H. Lipscomb, JGR 2001 203 !! 204 !! ** History : 205 !! largely inspired from CICE (c) W. H. Lipscomb and E.C. Hunke 206 !! 207 !! (01-2006) Martin Vancoppenolle, UCL-ASTR, translation from 208 !! CICE 209 !! (06-2006) Adaptation to include salt, age and types 210 !! (04-2007) Mass conservation checked 211 !!------------------------------------------------------------------ 212 !! * Arguments 213 214 INTEGER , INTENT (IN) :: & 215 klbnd , & ! Start thickness category index point 216 kubnd , & ! End point on which the the computation is applied 217 ntyp ! Number of the type used 218 219 !! * Local variables 220 INTEGER :: ji, & ! spatial dummy loop index 221 jj, & ! spatial dummy loop index 222 jl, & ! ice category dummy loop index 223 zji, zjj, & ! dummy indices used when changing coordinates 224 nd ! used for thickness categories 225 226 INTEGER , DIMENSION(jpi,jpj,jpl-1) :: & 227 zdonor ! donor category index 228 229 REAL(wp) :: & ! constant values 230 zeps = 1.0e-10 231 232 REAL(wp) :: & ! constant values for ice enthalpy 233 zindb , & 234 zareamin , & ! minimum tolerated area in a thickness category 235 zwk1, zwk2, & ! all the following are dummy arguments 236 zx1, zx2, zx3, & ! 237 zetamin , & ! minimum value of eta 238 zetamax , & ! maximum value of eta 239 zdh0 , & ! 240 zda0 , & ! 241 zdamax , & ! 242 zhimin 243 244 REAL(wp), DIMENSION(jpi,jpj,jpl) :: & 245 zdhice , & ! ice thickness increment 246 g0 , & ! coefficients for fitting the line of the ITD 247 g1 , & ! coefficients for fitting the line of the ITD 248 hL , & ! left boundary for the ITD for each thickness 249 hR , & ! left boundary for the ITD for each thickness 250 zht_i_o , & ! old ice thickness 251 dummy_es 252 253 REAL(wp), DIMENSION(jpi,jpj,jpl-1) :: & 254 zdaice , & ! local increment of ice area 255 zdvice ! local increment of ice volume 256 257 REAL(wp), DIMENSION(jpi,jpj,0:jpl) :: & 258 zhbnew ! new boundaries of ice categories 259 260 REAL(wp), DIMENSION(jpi,jpj) :: & 261 zhb0, zhb1 ! category boundaries for thinnes categories 262 263 REAL, DIMENSION(1:(jpi+1)*(jpj+1)) :: & 264 zvetamin, zvetamax ! maximum values for etas 265 266 INTEGER, DIMENSION(1:(jpi+1)*(jpj+1)) :: & 267 nind_i , & ! compressed indices for i/j directions 268 nind_j 269 270 INTEGER :: & 271 nbrem ! number of cells with ice to transfer 272 273 LOGICAL, DIMENSION(jpi,jpj) :: & !: 274 zremap_flag ! compute remapping or not ???? 275 276 REAL(wp) :: & ! constant values for ice enthalpy 277 zslope ! used to compute local thermodynamic "speeds" 278 279 REAL (wp), DIMENSION(jpi,jpj) :: & ! 280 vt_i_init, vt_i_final, & ! ice volume summed over categories 281 vt_s_init, vt_s_final, & ! snow volume summed over categories 282 et_i_init, et_i_final, & ! ice energy summed over categories 283 et_s_init, et_s_final ! snow energy summed over categories 284 285 CHARACTER (len = 15) :: fieldid 286 287 !!-- End of declarations 288 !!---------------------------------------------------------------------------------------------- 174 smv_i(:,:,:) = old_smv_i (:,:,:) 175 176 END SUBROUTINE lim_itd_th 177 ! 178 179 SUBROUTINE lim_itd_th_rem( klbnd, kubnd, ntyp, kt ) 180 !!------------------------------------------------------------------ 181 !! *** ROUTINE lim_itd_th_rem *** 182 !! ** Purpose : 183 !! This routine computes the redistribution of ice thickness 184 !! after thermodynamic growth of ice thickness 185 !! 186 !! ** Method : Linear remapping 187 !! 188 !! ** Arguments : 189 !! klbnd, kubnd : Starting and ending category index on which the 190 !! the computation is applied 191 !! 192 !! ** Inputs / Ouputs : (global commons) 193 !! 194 !! ** External : 195 !! 196 !! ** References : W.H. Lipscomb, JGR 2001 197 !! 198 !! ** History : 199 !! largely inspired from CICE (c) W. H. Lipscomb and E.C. Hunke 200 !! 201 !! (01-2006) Martin Vancoppenolle, UCL-ASTR, translation from 202 !! CICE 203 !! (06-2006) Adaptation to include salt, age and types 204 !! (04-2007) Mass conservation checked 205 !!------------------------------------------------------------------ 206 !! * Arguments 207 208 INTEGER , INTENT (IN) :: & 209 klbnd , & ! Start thickness category index point 210 kubnd , & ! End point on which the the computation is applied 211 ntyp , & ! Number of the type used 212 kt ! Ocean time step 213 214 !! * Local variables 215 INTEGER :: ji, & ! spatial dummy loop index 216 jj, & ! spatial dummy loop index 217 jl, & ! ice category dummy loop index 218 zji, zjj, & ! dummy indices used when changing coordinates 219 nd ! used for thickness categories 220 221 INTEGER , DIMENSION(jpi,jpj,jpl-1) :: & 222 zdonor ! donor category index 223 224 REAL(wp) :: & ! constant values 225 zeps = 1.0e-10 226 227 REAL(wp) :: & ! constant values for ice enthalpy 228 zindb , & 229 zareamin , & ! minimum tolerated area in a thickness category 230 zwk1, zwk2, & ! all the following are dummy arguments 231 zx1, zx2, zx3, & ! 232 zetamin , & ! minimum value of eta 233 zetamax , & ! maximum value of eta 234 zdh0 , & ! 235 zda0 , & ! 236 zdamax , & ! 237 zhimin 238 239 REAL(wp), DIMENSION(jpi,jpj,jpl) :: & 240 zdhice , & ! ice thickness increment 241 g0 , & ! coefficients for fitting the line of the ITD 242 g1 , & ! coefficients for fitting the line of the ITD 243 hL , & ! left boundary for the ITD for each thickness 244 hR , & ! left boundary for the ITD for each thickness 245 zht_i_o , & ! old ice thickness 246 dummy_es 247 248 REAL(wp), DIMENSION(jpi,jpj,jpl-1) :: & 249 zdaice , & ! local increment of ice area 250 zdvice ! local increment of ice volume 251 252 REAL(wp), DIMENSION(jpi,jpj,0:jpl) :: & 253 zhbnew ! new boundaries of ice categories 254 255 REAL(wp), DIMENSION(jpi,jpj) :: & 256 zhb0, zhb1 ! category boundaries for thinnes categories 257 258 REAL, DIMENSION(1:(jpi+1)*(jpj+1)) :: & 259 zvetamin, zvetamax ! maximum values for etas 260 261 INTEGER, DIMENSION(1:(jpi+1)*(jpj+1)) :: & 262 nind_i , & ! compressed indices for i/j directions 263 nind_j 264 265 INTEGER :: & 266 nbrem ! number of cells with ice to transfer 267 268 LOGICAL, DIMENSION(jpi,jpj) :: & !: 269 zremap_flag ! compute remapping or not ???? 270 271 REAL(wp) :: & ! constant values for ice enthalpy 272 zslope ! used to compute local thermodynamic "speeds" 273 274 REAL (wp), DIMENSION(jpi,jpj) :: & ! 275 vt_i_init, vt_i_final, & ! ice volume summed over categories 276 vt_s_init, vt_s_final, & ! snow volume summed over categories 277 et_i_init, et_i_final, & ! ice energy summed over categories 278 et_s_init, et_s_final ! snow energy summed over categories 279 280 CHARACTER (len = 15) :: fieldid 281 282 !!-- End of declarations 283 !!---------------------------------------------------------------------------------------------- 289 284 zhimin = 0.1 !minimum ice thickness tolerated by the model 290 285 zareamin = zeps !minimum area in thickness categories tolerated by the conceptors of the model 291 286 292 !!----------------------------------------------------------------------------------------------293 !! 0) Conservation checkand changes in each ice category294 !!----------------------------------------------------------------------------------------------287 !!---------------------------------------------------------------------------------------------- 288 !! 0) Conservation checkand changes in each ice category 289 !!---------------------------------------------------------------------------------------------- 295 290 IF ( con_i ) THEN 296 291 CALL lim_column_sum (jpl, v_i, vt_i_init) … … 300 295 CALL lim_column_sum (jpl, dummy_es(:,:,:) , et_s_init) 301 296 ENDIF 302 303 !!----------------------------------------------------------------------------------------------304 !! 1) Compute thickness and changes in each ice category305 !!----------------------------------------------------------------------------------------------306 IF (lwp) THEN307 WRITE(numout,*)308 WRITE(numout,*) 'lim_itd_th_rem : Remapping the ice thickness distribution'309 WRITE(numout,*) '~~~~~~~~~~~~~~~'310 WRITE(numout,*) ' klbnd : ', klbnd311 WRITE(numout,*) ' kubnd : ', kubnd312 WRITE(numout,*) ' ntyp : ', ntyp313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 !-----------------------------------------------------------------------------------------------331 ! 2) Compute fractional ice area in each grid cell332 !-----------------------------------------------------------------------------------------------297 298 !!---------------------------------------------------------------------------------------------- 299 !! 1) Compute thickness and changes in each ice category 300 !!---------------------------------------------------------------------------------------------- 301 IF (kt == nit000 .AND. lwp) THEN 302 WRITE(numout,*) 303 WRITE(numout,*) 'lim_itd_th_rem : Remapping the ice thickness distribution' 304 WRITE(numout,*) '~~~~~~~~~~~~~~~' 305 WRITE(numout,*) ' klbnd : ', klbnd 306 WRITE(numout,*) ' kubnd : ', kubnd 307 WRITE(numout,*) ' ntyp : ', ntyp 308 ENDIF 309 310 zdhice(:,:,:) = 0.0 311 DO jl = klbnd, kubnd 312 DO jj = 1, jpj 313 DO ji = 1, jpi 314 zindb = 1.0-MAX(0.0,SIGN(1.0,-a_i(ji,jj,jl))) !0 if no ice and 1 if yes 315 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX(a_i(ji,jj,jl),zeps) * zindb 316 zindb = 1.0-MAX(0.0,SIGN(1.0,-old_a_i(ji,jj,jl))) !0 if no ice and 1 if yes 317 zht_i_o(ji,jj,jl) = old_v_i(ji,jj,jl) / MAX(old_a_i(ji,jj,jl),zeps) * zindb 318 IF (a_i(ji,jj,jl).gt.1e-6) THEN 319 zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_o(ji,jj,jl) 320 ENDIF 321 END DO 322 END DO 323 END DO 324 325 !----------------------------------------------------------------------------------------------- 326 ! 2) Compute fractional ice area in each grid cell 327 !----------------------------------------------------------------------------------------------- 333 328 at_i(:,:) = 0.0 334 329 DO jl = klbnd, kubnd … … 340 335 END DO 341 336 342 !-----------------------------------------------------------------------------------------------343 ! 3) Identify grid cells with ice344 !-----------------------------------------------------------------------------------------------337 !----------------------------------------------------------------------------------------------- 338 ! 3) Identify grid cells with ice 339 !----------------------------------------------------------------------------------------------- 345 340 nbrem = 0 346 341 DO jj = 1, jpj … … 357 352 END DO !jj 358 353 359 !-----------------------------------------------------------------------------------------------360 ! 4) Compute new category boundaries361 !-----------------------------------------------------------------------------------------------354 !----------------------------------------------------------------------------------------------- 355 ! 4) Compute new category boundaries 356 !----------------------------------------------------------------------------------------------- 362 357 !- 4.1 Compute category boundaries 363 358 ! Tricky trick see limitd_me.F90 … … 374 369 ! 375 370 IF ( ( zht_i_o(zji,zjj,jl) .GT.zeps ) .AND. & 376 371 ( zht_i_o(zji,zjj,jl+1).GT.zeps ) ) THEN 377 372 !interpolate between adjacent category growth rates 378 373 zslope = ( zdhice(zji,zjj,jl+1) - zdhice(zji,zjj,jl) ) / & 379 374 ( zht_i_o (zji,zjj,jl+1) - zht_i_o (zji,zjj,jl) ) 380 375 zhbnew(zji,zjj,jl) = hi_max(jl) + zdhice(zji,zjj,jl) + & 381 376 zslope * ( hi_max(jl) - zht_i_o(zji,zjj,jl) ) 382 377 ELSEIF (zht_i_o(zji,zjj,jl).gt.zeps) THEN 383 378 zhbnew(zji,zjj,jl) = hi_max(jl) + zdhice(zji,zjj,jl) … … 391 386 ! jl 392 387 393 !- 4.2 Check that each zhbnew lies between adjacent values of ice thickness388 !- 4.2 Check that each zhbnew lies between adjacent values of ice thickness 394 389 DO ji = 1, nbrem 395 390 ! jl, ji … … 398 393 ! jl, ji 399 394 IF ( ( a_i(zji,zjj,jl) .GT.zeps) .AND. & 400 395 ( ht_i(zji,zjj,jl).GE. zhbnew(zji,zjj,jl) ) & 401 396 ) THEN 402 397 zremap_flag(zji,zjj) = .false. 403 398 ELSEIF ( ( a_i(zji,zjj,jl+1) .GT. zeps ) .AND. & 404 405 399 ( ht_i(zji,zjj,jl+1).LE. zhbnew(zji,zjj,jl) ) & 400 ) THEN 406 401 zremap_flag(zji,zjj) = .false. 407 402 ENDIF 408 403 409 !- 4.3 Check that each zhbnew does not exceed maximal values hi_max404 !- 4.3 Check that each zhbnew does not exceed maximal values hi_max 410 405 ! jl, ji 411 406 IF (zhbnew(zji,zjj,jl).gt.hi_max(jl+1)) THEN … … 420 415 ! ji 421 416 END DO !jl 422 423 !-----------------------------------------------------------------------------------------------424 ! 5) Identify cells where ITD is to be remapped425 !-----------------------------------------------------------------------------------------------426 nbrem = 0427 DO jj = 1, jpj428 DO ji = 1, jpi429 IF ( zremap_flag(ji,jj) ) THEN430 nbrem = nbrem + 1431 nind_i(nbrem) = ji432 nind_j(nbrem) = jj433 ENDIF434 END DO !ji435 END DO !jj436 437 !-----------------------------------------------------------------------------------------------438 ! 6) Fill arrays with lowermost / uppermost boundaries of 'new' categories439 !-----------------------------------------------------------------------------------------------440 DO jj = 1, jpj441 DO ji = 1, jpi442 zhb0(ji,jj) = hi_max_typ(0,ntyp) ! 0eme443 zhb1(ji,jj) = hi_max_typ(1,ntyp) ! 1er444 445 zhbnew(ji,jj,klbnd-1) = 0.0446 447 IF ( a_i(ji,jj,kubnd) .GT. zeps ) THEN448 zhbnew(ji,jj,kubnd) = 3.0*ht_i(ji,jj,kubnd) - 2.0*zhbnew(ji,jj,kubnd-1)449 ELSE450 zhbnew(ji,jj,kubnd) = hi_max(kubnd)451 ENDIF452 453 IF ( zhbnew(ji,jj,kubnd) .LT. hi_max(kubnd-1) ) &454 zhbnew(ji,jj,kubnd) = hi_max(kubnd-1)455 456 END DO !jj457 END DO !jj458 459 !-----------------------------------------------------------------------------------------------460 ! 7) Compute g(h)461 !-----------------------------------------------------------------------------------------------462 !- 7.1 g(h) for category 1 at start of time step463 CALL lim_itd_fitline(klbnd, zhb0, zhb1, zht_i_o(:,:,klbnd), &464 465 466 467 !- 7.2 Area lost due to melting of thin ice (first category, klbnd)468 DO ji = 1, nbrem469 zji = nind_i(ji)470 zjj = nind_j(ji)471 472 !ji473 IF (a_i(zji,zjj,klbnd) .gt. zeps) THEN474 zdh0 = zdhice(zji,zjj,klbnd) !decrease of ice thickness in the lower category475 ! ji, a_i > zeps476 IF (zdh0 .lt. 0.0) THEN !remove area from category 1477 ! ji, a_i > zeps; zdh0 < 0478 zdh0 = MIN(-zdh0,hi_max(klbnd))479 480 !Integrate g(1) from 0 to dh0 to estimate area melted481 zetamax = MIN(zdh0,hR(zji,zjj,klbnd)) - hL(zji,zjj,klbnd)482 IF (zetamax.gt.0.0) THEN483 zx1 = zetamax484 zx2 = 0.5 * zetamax*zetamax485 zda0 = g1(zji,zjj,klbnd) * zx2 + g0(zji,zjj,klbnd) * zx1 !ice area removed486 ! Constrain new thickness <= ht_i487 zdamax = a_i(zji,zjj,klbnd) * &488 489 !ice area lost due to melting of thin ice490 zda0 = MIN(zda0, zdamax)491 492 ! Remove area, conserving volume493 ht_i(zji,zjj,klbnd) = ht_i(zji,zjj,klbnd) &494 495 a_i(zji,zjj,klbnd) = a_i(zji,zjj,klbnd) - zda0496 v_i(zji,zjj,klbnd) = a_i(zji,zjj,klbnd)*ht_i(zji,zjj,klbnd)497 ENDIF ! zetamax > 0498 ! ji, a_i > zeps499 500 ELSE ! if ice accretion501 ! ji, a_i > zeps; zdh0 > 0502 IF ( ntyp .EQ. 1 ) zhbnew(zji,zjj,klbnd-1) = MIN(zdh0,hi_max(klbnd))503 ! zhbnew was 0, and is shifted to the right to account for thin ice504 ! growth in openwater (F0 = f1)505 IF ( ntyp .NE. 1 ) zhbnew(zji,zjj,0) = 0506 ! in other types there is507 ! no open water growth (F0 = 0)508 ENDIF ! zdh0509 510 ! a_i > zeps511 ENDIF ! a_i > zeps512 513 END DO ! ji514 515 !- 7.3 g(h) for each thickness category516 DO jl = klbnd, kubnd517 CALL lim_itd_fitline(jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl), &518 519 520 END DO521 522 !-----------------------------------------------------------------------------------------------523 ! 8) Compute area and volume to be shifted across each boundary524 !-----------------------------------------------------------------------------------------------525 526 DO jl = klbnd, kubnd - 1527 DO jj = 1, jpj528 DO ji = 1, jpi529 zdonor(ji,jj,jl) = 0530 zdaice(ji,jj,jl) = 0.0531 zdvice(ji,jj,jl) = 0.0532 END DO533 END DO534 535 DO ji = 1, nbrem536 zji = nind_i(ji)537 zjj = nind_j(ji)538 539 IF (zhbnew(zji,zjj,jl) .gt. hi_max(jl)) THEN ! transfer from jl to jl+1540 541 ! left and right integration limits in eta space542 zvetamin(ji) = MAX(hi_max(jl), hL(zji,zjj,jl)) - hL(zji,zjj,jl)543 zvetamax(ji) = MIN(zhbnew(zji,zjj,jl), hR(zji,zjj,jl)) - hL(zji,zjj,jl)544 zdonor(zji,zjj,jl) = jl545 546 ELSE ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl547 548 ! left and right integration limits in eta space549 zvetamin(ji) = 0.0550 zvetamax(ji) = MIN(hi_max(jl), hR(zji,zjj,jl+1)) - hL(zji,zjj,jl+1)551 zdonor(zji,zjj,jl) = jl + 1552 553 ENDIF ! zhbnew(jl) > hi_max(jl)554 555 zetamax = MAX(zvetamax(ji), zvetamin(ji)) ! no transfer if etamax < etamin556 zetamin = zvetamin(ji)557 558 zx1 = zetamax - zetamin559 zwk1 = zetamin*zetamin560 zwk2 = zetamax*zetamax561 zx2 = 0.5 * (zwk2 - zwk1)562 zwk1 = zwk1 * zetamin563 zwk2 = zwk2 * zetamax564 zx3 = 1.0/3.0 * (zwk2 - zwk1)565 nd = zdonor(zji,zjj,jl)566 zdaice(zji,zjj,jl) = g1(zji,zjj,nd)*zx2 + g0(zji,zjj,nd)*zx1567 zdvice(zji,zjj,jl) = g1(zji,zjj,nd)*zx3 + g0(zji,zjj,nd)*zx2 + &568 569 570 END DO ! ji571 END DO ! jl klbnd -> kubnd - 1572 573 !!----------------------------------------------------------------------------------------------574 !! 9) Shift ice between categories575 !!----------------------------------------------------------------------------------------------576 CALL lim_itd_shiftice ( klbnd, kubnd, zdonor, zdaice, zdvice )577 578 !!----------------------------------------------------------------------------------------------579 !! 10) Make sure ht_i >= minimum ice thickness hi_min580 !!----------------------------------------------------------------------------------------------581 582 DO ji = 1, nbrem583 zji = nind_i(ji)584 zjj = nind_j(ji)585 IF ( ( zhimin .GT. 0.0 ) .AND. &586 587 ) THEN588 a_i(zji,zjj,1) = a_i(zji,zjj,1) * ht_i(zji,zjj,1) / zhimin589 ht_i(zji,zjj,1) = zhimin590 v_i(zji,zjj,1) = a_i(zji,zjj,1)*ht_i(zji,zjj,1)591 ENDIF592 END DO !ji593 594 !!----------------------------------------------------------------------------------------------595 !! 11) Conservation check596 !!----------------------------------------------------------------------------------------------417 418 !----------------------------------------------------------------------------------------------- 419 ! 5) Identify cells where ITD is to be remapped 420 !----------------------------------------------------------------------------------------------- 421 nbrem = 0 422 DO jj = 1, jpj 423 DO ji = 1, jpi 424 IF ( zremap_flag(ji,jj) ) THEN 425 nbrem = nbrem + 1 426 nind_i(nbrem) = ji 427 nind_j(nbrem) = jj 428 ENDIF 429 END DO !ji 430 END DO !jj 431 432 !----------------------------------------------------------------------------------------------- 433 ! 6) Fill arrays with lowermost / uppermost boundaries of 'new' categories 434 !----------------------------------------------------------------------------------------------- 435 DO jj = 1, jpj 436 DO ji = 1, jpi 437 zhb0(ji,jj) = hi_max_typ(0,ntyp) ! 0eme 438 zhb1(ji,jj) = hi_max_typ(1,ntyp) ! 1er 439 440 zhbnew(ji,jj,klbnd-1) = 0.0 441 442 IF ( a_i(ji,jj,kubnd) .GT. zeps ) THEN 443 zhbnew(ji,jj,kubnd) = 3.0*ht_i(ji,jj,kubnd) - 2.0*zhbnew(ji,jj,kubnd-1) 444 ELSE 445 zhbnew(ji,jj,kubnd) = hi_max(kubnd) 446 ENDIF 447 448 IF ( zhbnew(ji,jj,kubnd) .LT. hi_max(kubnd-1) ) & 449 zhbnew(ji,jj,kubnd) = hi_max(kubnd-1) 450 451 END DO !jj 452 END DO !jj 453 454 !----------------------------------------------------------------------------------------------- 455 ! 7) Compute g(h) 456 !----------------------------------------------------------------------------------------------- 457 !- 7.1 g(h) for category 1 at start of time step 458 CALL lim_itd_fitline(klbnd, zhb0, zhb1, zht_i_o(:,:,klbnd), & 459 g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd), & 460 hR(:,:,klbnd), zremap_flag) 461 462 !- 7.2 Area lost due to melting of thin ice (first category, klbnd) 463 DO ji = 1, nbrem 464 zji = nind_i(ji) 465 zjj = nind_j(ji) 466 467 !ji 468 IF (a_i(zji,zjj,klbnd) .gt. zeps) THEN 469 zdh0 = zdhice(zji,zjj,klbnd) !decrease of ice thickness in the lower category 470 ! ji, a_i > zeps 471 IF (zdh0 .lt. 0.0) THEN !remove area from category 1 472 ! ji, a_i > zeps; zdh0 < 0 473 zdh0 = MIN(-zdh0,hi_max(klbnd)) 474 475 !Integrate g(1) from 0 to dh0 to estimate area melted 476 zetamax = MIN(zdh0,hR(zji,zjj,klbnd)) - hL(zji,zjj,klbnd) 477 IF (zetamax.gt.0.0) THEN 478 zx1 = zetamax 479 zx2 = 0.5 * zetamax*zetamax 480 zda0 = g1(zji,zjj,klbnd) * zx2 + g0(zji,zjj,klbnd) * zx1 !ice area removed 481 ! Constrain new thickness <= ht_i 482 zdamax = a_i(zji,zjj,klbnd) * & 483 (1.0 - ht_i(zji,zjj,klbnd)/zht_i_o(zji,zjj,klbnd)) ! zdamax > 0 484 !ice area lost due to melting of thin ice 485 zda0 = MIN(zda0, zdamax) 486 487 ! Remove area, conserving volume 488 ht_i(zji,zjj,klbnd) = ht_i(zji,zjj,klbnd) & 489 * a_i(zji,zjj,klbnd) / ( a_i(zji,zjj,klbnd) - zda0 ) 490 a_i(zji,zjj,klbnd) = a_i(zji,zjj,klbnd) - zda0 491 v_i(zji,zjj,klbnd) = a_i(zji,zjj,klbnd)*ht_i(zji,zjj,klbnd) 492 ENDIF ! zetamax > 0 493 ! ji, a_i > zeps 494 495 ELSE ! if ice accretion 496 ! ji, a_i > zeps; zdh0 > 0 497 IF ( ntyp .EQ. 1 ) zhbnew(zji,zjj,klbnd-1) = MIN(zdh0,hi_max(klbnd)) 498 ! zhbnew was 0, and is shifted to the right to account for thin ice 499 ! growth in openwater (F0 = f1) 500 IF ( ntyp .NE. 1 ) zhbnew(zji,zjj,0) = 0 501 ! in other types there is 502 ! no open water growth (F0 = 0) 503 ENDIF ! zdh0 504 505 ! a_i > zeps 506 ENDIF ! a_i > zeps 507 508 END DO ! ji 509 510 !- 7.3 g(h) for each thickness category 511 DO jl = klbnd, kubnd 512 CALL lim_itd_fitline(jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl), & 513 g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), & 514 zremap_flag) 515 END DO 516 517 !----------------------------------------------------------------------------------------------- 518 ! 8) Compute area and volume to be shifted across each boundary 519 !----------------------------------------------------------------------------------------------- 520 521 DO jl = klbnd, kubnd - 1 522 DO jj = 1, jpj 523 DO ji = 1, jpi 524 zdonor(ji,jj,jl) = 0 525 zdaice(ji,jj,jl) = 0.0 526 zdvice(ji,jj,jl) = 0.0 527 END DO 528 END DO 529 530 DO ji = 1, nbrem 531 zji = nind_i(ji) 532 zjj = nind_j(ji) 533 534 IF (zhbnew(zji,zjj,jl) .gt. hi_max(jl)) THEN ! transfer from jl to jl+1 535 536 ! left and right integration limits in eta space 537 zvetamin(ji) = MAX(hi_max(jl), hL(zji,zjj,jl)) - hL(zji,zjj,jl) 538 zvetamax(ji) = MIN(zhbnew(zji,zjj,jl), hR(zji,zjj,jl)) - hL(zji,zjj,jl) 539 zdonor(zji,zjj,jl) = jl 540 541 ELSE ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl 542 543 ! left and right integration limits in eta space 544 zvetamin(ji) = 0.0 545 zvetamax(ji) = MIN(hi_max(jl), hR(zji,zjj,jl+1)) - hL(zji,zjj,jl+1) 546 zdonor(zji,zjj,jl) = jl + 1 547 548 ENDIF ! zhbnew(jl) > hi_max(jl) 549 550 zetamax = MAX(zvetamax(ji), zvetamin(ji)) ! no transfer if etamax < etamin 551 zetamin = zvetamin(ji) 552 553 zx1 = zetamax - zetamin 554 zwk1 = zetamin*zetamin 555 zwk2 = zetamax*zetamax 556 zx2 = 0.5 * (zwk2 - zwk1) 557 zwk1 = zwk1 * zetamin 558 zwk2 = zwk2 * zetamax 559 zx3 = 1.0/3.0 * (zwk2 - zwk1) 560 nd = zdonor(zji,zjj,jl) 561 zdaice(zji,zjj,jl) = g1(zji,zjj,nd)*zx2 + g0(zji,zjj,nd)*zx1 562 zdvice(zji,zjj,jl) = g1(zji,zjj,nd)*zx3 + g0(zji,zjj,nd)*zx2 + & 563 zdaice(zji,zjj,jl)*hL(zji,zjj,nd) 564 565 END DO ! ji 566 END DO ! jl klbnd -> kubnd - 1 567 568 !!---------------------------------------------------------------------------------------------- 569 !! 9) Shift ice between categories 570 !!---------------------------------------------------------------------------------------------- 571 CALL lim_itd_shiftice ( klbnd, kubnd, zdonor, zdaice, zdvice ) 572 573 !!---------------------------------------------------------------------------------------------- 574 !! 10) Make sure ht_i >= minimum ice thickness hi_min 575 !!---------------------------------------------------------------------------------------------- 576 577 DO ji = 1, nbrem 578 zji = nind_i(ji) 579 zjj = nind_j(ji) 580 IF ( ( zhimin .GT. 0.0 ) .AND. & 581 ( ( a_i(zji,zjj,1) .GT. zeps ) .AND. ( ht_i(zji,zjj,1) .LT. zhimin ) ) & 582 ) THEN 583 a_i(zji,zjj,1) = a_i(zji,zjj,1) * ht_i(zji,zjj,1) / zhimin 584 ht_i(zji,zjj,1) = zhimin 585 v_i(zji,zjj,1) = a_i(zji,zjj,1)*ht_i(zji,zjj,1) 586 ENDIF 587 END DO !ji 588 589 !!---------------------------------------------------------------------------------------------- 590 !! 11) Conservation check 591 !!---------------------------------------------------------------------------------------------- 597 592 IF ( con_i ) THEN 598 593 CALL lim_column_sum (jpl, v_i, vt_i_final) … … 614 609 ENDIF 615 610 616 END SUBROUTINE lim_itd_th_rem 617 618 !!---------------------------------------------------------------------------------------------- 619 !!---------------------------------------------------------------------------------------------- 620 621 SUBROUTINE lim_itd_fitline(num_cat, HbL, Hbr, hice, g0, g1, hL, hR, zremap_flag ) 622 623 !!------------------------------------------------------------------ 624 !! *** ROUTINE lim_itd_fitline *** 625 !! ** Purpose : 626 !! fit g(h) with a line using area, volume constraints 627 !! 628 !! ** Method : 629 !! Fit g(h) with a line, satisfying area and volume constraints. 630 !! To reduce roundoff errors caused by large values of g0 and g1, 631 !! we actually compute g(eta), where eta = h - hL, and hL is the 632 !! left boundary. 633 !! 634 !! ** Arguments : 635 !! 636 !! ** Inputs / Ouputs : (global commons) 637 !! 638 !! ** External : 639 !! 640 !! ** References : 641 !! 642 !! ** History : 643 !! authors: William H. Lipscomb, LANL, Elizabeth C. Hunke, LANL 644 !! (01-2006) Martin Vancoppenolle 645 !! 646 !!------------------------------------------------------------------ 647 !! * Arguments 611 END SUBROUTINE lim_itd_th_rem 612 ! 613 614 SUBROUTINE lim_itd_fitline(num_cat, HbL, Hbr, hice, g0, g1, hL, hR, zremap_flag ) 615 616 !!------------------------------------------------------------------ 617 !! *** ROUTINE lim_itd_fitline *** 618 !! ** Purpose : 619 !! fit g(h) with a line using area, volume constraints 620 !! 621 !! ** Method : 622 !! Fit g(h) with a line, satisfying area and volume constraints. 623 !! To reduce roundoff errors caused by large values of g0 and g1, 624 !! we actually compute g(eta), where eta = h - hL, and hL is the 625 !! left boundary. 626 !! 627 !! ** Arguments : 628 !! 629 !! ** Inputs / Ouputs : (global commons) 630 !! 631 !! ** External : 632 !! 633 !! ** References : 634 !! 635 !! ** History : 636 !! authors: William H. Lipscomb, LANL, Elizabeth C. Hunke, LANL 637 !! (01-2006) Martin Vancoppenolle 638 !! 639 !!------------------------------------------------------------------ 640 !! * Arguments 648 641 649 642 INTEGER, INTENT(in) :: num_cat ! category index … … 674 667 675 668 REAL(wp) :: & ! constant values 676 669 zeps = 1.0e-10 677 670 678 671 zacrith = 1.0e-6 679 !!-- End of declarations680 !!----------------------------------------------------------------------------------------------672 !!-- End of declarations 673 !!---------------------------------------------------------------------------------------------- 681 674 682 675 DO jj = 1, jpj … … 684 677 685 678 IF ( zremap_flag(ji,jj) .AND. a_i(ji,jj,num_cat) .gt. zacrith & 686 687 688 ! Initialize hL and hR679 .AND. hice(ji,jj) .GT. 0.0 ) THEN 680 681 ! Initialize hL and hR 689 682 690 683 hL(ji,jj) = HbL(ji,jj) 691 684 hR(ji,jj) = HbR(ji,jj) 692 685 693 ! Change hL or hR if hice falls outside central third of range686 ! Change hL or hR if hice falls outside central third of range 694 687 695 688 zh13 = 1.0/3.0 * (2.0*hL(ji,jj) + hR(ji,jj)) … … 702 695 ENDIF 703 696 704 ! Compute coefficients of g(eta) = g0 + g1*eta705 697 ! Compute coefficients of g(eta) = g0 + g1*eta 698 706 699 zdhr = 1.0 / (hR(ji,jj) - hL(ji,jj)) 707 700 zwk1 = 6.0 * a_i(ji,jj,num_cat) * zdhr … … 722 715 END DO ! jj 723 716 724 END SUBROUTINE lim_itd_fitline 725 726 !---------------------------------------------------------------------------------------------- 727 !---------------------------------------------------------------------------------------------- 728 729 SUBROUTINE lim_itd_shiftice (klbnd, kubnd, zdonor, zdaice, zdvice) 730 !!------------------------------------------------------------------ 731 !! *** ROUTINE lim_itd_shiftice *** 732 !! ** Purpose : shift ice across category boundaries, conserving everything 733 !! ( area, volume, energy, age*vol, and mass of salt ) 734 !! 735 !! ** Method : 736 !! 737 !! ** Arguments : 738 !! 739 !! ** Inputs / Ouputs : (global commons) 740 !! 741 !! ** External : 742 !! 743 !! ** References : 744 !! 745 !! ** History : 746 !! authors: William H. Lipscomb, LANL, Elizabeth C. Hunke, LANL 747 !! (01-2006) Martin Vancoppenolle 748 !! 749 !!------------------------------------------------------------------ 750 !! * Arguments 717 END SUBROUTINE lim_itd_fitline 718 ! 719 720 SUBROUTINE lim_itd_shiftice (klbnd, kubnd, zdonor, zdaice, zdvice) 721 !!------------------------------------------------------------------ 722 !! *** ROUTINE lim_itd_shiftice *** 723 !! ** Purpose : shift ice across category boundaries, conserving everything 724 !! ( area, volume, energy, age*vol, and mass of salt ) 725 !! 726 !! ** Method : 727 !! 728 !! ** Arguments : 729 !! 730 !! ** Inputs / Ouputs : (global commons) 731 !! 732 !! ** External : 733 !! 734 !! ** References : 735 !! 736 !! ** History : 737 !! authors: William H. Lipscomb, LANL, Elizabeth C. Hunke, LANL 738 !! (01-2006) Martin Vancoppenolle 739 !! 740 !!------------------------------------------------------------------ 741 !! * Arguments 751 742 752 743 INTEGER , INTENT (IN) :: & 753 754 744 klbnd , & ! Start thickness category index point 745 kubnd ! End point on which the the computation is applied 755 746 756 747 INTEGER , DIMENSION(jpi,jpj,jpl-1), INTENT(IN) :: & … … 792 783 793 784 LOGICAL :: & 794 zdaice_negative , & ! true if daice < -puny795 zdvice_negative , & ! true if dvice < -puny796 zdaice_greater_aicen , & ! true if daice > aicen797 zdvice_greater_vicen ! true if dvice > vicen798 799 800 801 802 !!-- End of declarations803 804 !----------------------------------------------------------------------------------------------805 ! 1) Define a variable equal to a_i*T_su806 !----------------------------------------------------------------------------------------------785 zdaice_negative , & ! true if daice < -puny 786 zdvice_negative , & ! true if dvice < -puny 787 zdaice_greater_aicen , & ! true if daice > aicen 788 zdvice_greater_vicen ! true if dvice > vicen 789 790 REAL(wp) :: & ! constant values 791 zeps = 1.0e-10 792 793 !!-- End of declarations 794 795 !---------------------------------------------------------------------------------------------- 796 ! 1) Define a variable equal to a_i*T_su 797 !---------------------------------------------------------------------------------------------- 807 798 808 799 DO jl = klbnd, kubnd … … 814 805 END DO ! jl 815 806 816 !----------------------------------------------------------------------------------------------817 ! 2) Check for daice or dvice out of range, allowing for roundoff error818 !----------------------------------------------------------------------------------------------807 !---------------------------------------------------------------------------------------------- 808 ! 2) Check for daice or dvice out of range, allowing for roundoff error 809 !---------------------------------------------------------------------------------------------- 819 810 ! Note: zdaice < 0 or zdvice < 0 usually happens when category jl 820 811 ! has a small area, with h(n) very close to a boundary. Then … … 834 825 DO ji = 1, jpi 835 826 836 IF (zdonor(ji,jj,jl) .GT. 0) THEN 837 jl1 = zdonor(ji,jj,jl) 838 839 IF (zdaice(ji,jj,jl) .LT. 0.0) THEN 840 IF (zdaice(ji,jj,jl) .GT. -zeps) THEN 841 IF ( ( jl1.EQ.jl .AND. ht_i(ji,jj,jl1) .GT. hi_max(jl) ) & 842 .OR. & 843 ( jl1.EQ.jl+1 .AND. ht_i(ji,jj,jl1) .LE. hi_max(jl) ) & 844 ) THEN 845 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category 846 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 827 IF (zdonor(ji,jj,jl) .GT. 0) THEN 828 jl1 = zdonor(ji,jj,jl) 829 830 IF (zdaice(ji,jj,jl) .LT. 0.0) THEN 831 IF (zdaice(ji,jj,jl) .GT. -zeps) THEN 832 IF ( ( jl1.EQ.jl .AND. ht_i(ji,jj,jl1) .GT. hi_max(jl) ) & 833 .OR. & 834 ( jl1.EQ.jl+1 .AND. ht_i(ji,jj,jl1) .LE. hi_max(jl) ) & 835 ) THEN 836 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category 837 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 838 ELSE 839 zdaice(ji,jj,jl) = 0.0 ! shift no ice 840 zdvice(ji,jj,jl) = 0.0 841 ENDIF 847 842 ELSE 848 zdaice(ji,jj,jl) = 0.0 ! shift no ice 849 zdvice(ji,jj,jl) = 0.0 843 zdaice_negative = .true. 850 844 ENDIF 851 ELSE852 zdaice_negative = .true.853 845 ENDIF 854 ENDIF 855 856 IF (zdvice(ji,jj,jl) .LT. 0.0) THEN 857 IF (zdvice(ji,jj,jl) .GT. -zeps ) THEN 858 IF ( ( jl1.EQ.jl .AND. ht_i(ji,jj,jl1).GT.hi_max(jl) ) & 859 .OR. & 860 ( jl1.EQ.jl+1 .AND. ht_i(ji,jj,jl1) .LE. hi_max(jl) ) & 861 ) THEN 862 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category 846 847 IF (zdvice(ji,jj,jl) .LT. 0.0) THEN 848 IF (zdvice(ji,jj,jl) .GT. -zeps ) THEN 849 IF ( ( jl1.EQ.jl .AND. ht_i(ji,jj,jl1).GT.hi_max(jl) ) & 850 .OR. & 851 ( jl1.EQ.jl+1 .AND. ht_i(ji,jj,jl1) .LE. hi_max(jl) ) & 852 ) THEN 853 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category 854 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 855 ELSE 856 zdaice(ji,jj,jl) = 0.0 ! shift no ice 857 zdvice(ji,jj,jl) = 0.0 858 ENDIF 859 ELSE 860 zdvice_negative = .true. 861 ENDIF 862 ENDIF 863 864 ! If daice is close to aicen, set daice = aicen. 865 IF (zdaice(ji,jj,jl) .GT. a_i(ji,jj,jl1) - zeps ) THEN 866 IF (zdaice(ji,jj,jl) .LT. a_i(ji,jj,jl1)+zeps) THEN 867 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 863 868 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 864 869 ELSE 865 zdaice(ji,jj,jl) = 0.0 ! shift no ice 866 zdvice(ji,jj,jl) = 0.0 870 zdaice_greater_aicen = .true. 867 871 ENDIF 868 ELSE869 zdvice_negative = .true.870 872 ENDIF 871 ENDIF 872 873 ! If daice is close to aicen, set daice = aicen. 874 IF (zdaice(ji,jj,jl) .GT. a_i(ji,jj,jl1) - zeps ) THEN 875 IF (zdaice(ji,jj,jl) .LT. a_i(ji,jj,jl1)+zeps) THEN 876 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 877 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 878 ELSE 879 zdaice_greater_aicen = .true. 873 874 IF (zdvice(ji,jj,jl) .GT. v_i(ji,jj,jl1)-zeps) THEN 875 IF (zdvice(ji,jj,jl) .LT. v_i(ji,jj,jl1)+zeps) THEN 876 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 877 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 878 ELSE 879 zdvice_greater_vicen = .true. 880 ENDIF 880 881 ENDIF 881 ENDIF 882 883 IF (zdvice(ji,jj,jl) .GT. v_i(ji,jj,jl1)-zeps) THEN 884 IF (zdvice(ji,jj,jl) .LT. v_i(ji,jj,jl1)+zeps) THEN 885 zdaice(ji,jj,jl) = a_i(ji,jj,jl1) 886 zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 887 ELSE 888 zdvice_greater_vicen = .true. 889 ENDIF 890 ENDIF 891 892 ENDIF ! donor > 0 893 END DO ! i 882 883 ENDIF ! donor > 0 884 END DO ! i 894 885 END DO ! j 895 886 896 887 END DO !jl 897 888 898 !-------------------------------------------------------------------------------899 ! 3) Transfer volume and energy between categories900 !-------------------------------------------------------------------------------889 !------------------------------------------------------------------------------- 890 ! 3) Transfer volume and energy between categories 891 !------------------------------------------------------------------------------- 901 892 902 893 DO jl = klbnd, kubnd - 1 … … 1012 1003 DO jl = klbnd, kubnd 1013 1004 DO jj = 1, jpj 1014 DO ji = 1, jpi1015 IF ( a_i(ji,jj,jl) .GT. zeps ) THEN1016 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl)1017 t_su(ji,jj,jl) = zaTsfn(ji,jj,jl) / a_i(ji,jj,jl)1018 zindsn = 1.0 - MAX(0.0,SIGN(1.0,-v_s(ji,jj,jl))) !0 if no ice and 1 if yes1019 ELSE1020 ht_i(ji,jj,jl) = 0.01021 t_su(ji,jj,jl) = rtt1022 ENDIF1023 END DO ! ji1005 DO ji = 1, jpi 1006 IF ( a_i(ji,jj,jl) .GT. zeps ) THEN 1007 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 1008 t_su(ji,jj,jl) = zaTsfn(ji,jj,jl) / a_i(ji,jj,jl) 1009 zindsn = 1.0 - MAX(0.0,SIGN(1.0,-v_s(ji,jj,jl))) !0 if no ice and 1 if yes 1010 ELSE 1011 ht_i(ji,jj,jl) = 0.0 1012 t_su(ji,jj,jl) = rtt 1013 ENDIF 1014 END DO ! ji 1024 1015 END DO ! jj 1025 1016 END DO ! jl 1026 1017 1027 END SUBROUTINE lim_itd_shiftice 1028 1029 !---------------------------------------------------------------------------------------- 1030 !---------------------------------------------------------------------------------------- 1031 1032 SUBROUTINE lim_itd_th_reb(klbnd, kubnd, ntyp) 1033 !!------------------------------------------------------------------ 1034 !! *** ROUTINE lim_itd_th_reb *** 1035 !! ** Purpose : rebin - rebins thicknesses into defined categories 1036 !! 1037 !! ** Method : 1038 !! 1039 !! ** Arguments : 1040 !! 1041 !! ** Inputs / Ouputs : (global commons) 1042 !! 1043 !! ** External : 1044 !! 1045 !! ** References : 1046 !! 1047 !! ** History : (2005) Translation from CICE 1048 !! (2006) Adaptation to include salt, age and types 1049 !! (2007) Mass conservation checked 1050 !! 1051 !! authors: William H. Lipscomb, LANL, Elizabeth C. Hunke, LANL 1052 !! (01-2006) Martin Vancoppenolle (adaptation) 1053 !! 1054 !!------------------------------------------------------------------ 1055 !! * Arguments 1018 END SUBROUTINE lim_itd_shiftice 1019 ! 1020 1021 SUBROUTINE lim_itd_th_reb(klbnd, kubnd, ntyp) 1022 !!------------------------------------------------------------------ 1023 !! *** ROUTINE lim_itd_th_reb *** 1024 !! ** Purpose : rebin - rebins thicknesses into defined categories 1025 !! 1026 !! ** Method : 1027 !! 1028 !! ** Arguments : 1029 !! 1030 !! ** Inputs / Ouputs : (global commons) 1031 !! 1032 !! ** External : 1033 !! 1034 !! ** References : 1035 !! 1036 !! ** History : (2005) Translation from CICE 1037 !! (2006) Adaptation to include salt, age and types 1038 !! (2007) Mass conservation checked 1039 !! 1040 !! authors: William H. Lipscomb, LANL, Elizabeth C. Hunke, LANL 1041 !! (01-2006) Martin Vancoppenolle (adaptation) 1042 !! 1043 !!------------------------------------------------------------------ 1044 !! * Arguments 1056 1045 INTEGER , INTENT (in) :: & 1057 1058 1059 1046 klbnd , & ! Start thickness category index point 1047 kubnd , & ! End point on which the the computation is applied 1048 ntyp ! number of the ice type involved in the rebinning process 1060 1049 1061 1050 INTEGER :: & … … 1081 1070 vt_s_init, vt_s_final ! snow volume summed over categories 1082 1071 1083 1084 1085 !!-- End of declarations1086 !------------------------------------------------------------------------------1087 1088 ! ! conservation check1072 CHARACTER (len = 15) :: fieldid 1073 1074 !!-- End of declarations 1075 !------------------------------------------------------------------------------ 1076 1077 ! ! conservation check 1089 1078 IF ( con_i ) THEN 1090 1079 CALL lim_column_sum (jpl, v_i, vt_i_init) … … 1092 1081 ENDIF 1093 1082 1094 !1095 !------------------------------------------------------------------------------1096 ! 1) Compute ice thickness.1097 !------------------------------------------------------------------------------1083 ! 1084 !------------------------------------------------------------------------------ 1085 ! 1) Compute ice thickness. 1086 !------------------------------------------------------------------------------ 1098 1087 DO jl = klbnd, kubnd 1099 1088 DO jj = 1, jpj 1100 DO ji = 1, jpi1101 IF (a_i(ji,jj,jl) .GT. zeps) THEN1102 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl)1103 ELSE1104 ht_i(ji,jj,jl) = 0.01105 ENDIF1106 END DO ! i1089 DO ji = 1, jpi 1090 IF (a_i(ji,jj,jl) .GT. zeps) THEN 1091 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 1092 ELSE 1093 ht_i(ji,jj,jl) = 0.0 1094 ENDIF 1095 END DO ! i 1107 1096 END DO ! j 1108 1097 END DO ! n 1109 1098 1110 !------------------------------------------------------------------------------1111 ! 2) Make sure thickness of cat klbnd is at least hi_max_typ(klbnd)1112 !------------------------------------------------------------------------------1099 !------------------------------------------------------------------------------ 1100 ! 2) Make sure thickness of cat klbnd is at least hi_max_typ(klbnd) 1101 !------------------------------------------------------------------------------ 1113 1102 DO jj = 1, jpj 1114 DO ji = 1, jpi 1115 1116 IF (a_i(ji,jj,klbnd) > zeps) THEN 1117 IF (ht_i(ji,jj,klbnd) .LE. hi_max_typ(0,ntyp) .AND. hi_max_typ(0,ntyp) .GT. 0.0 ) THEN 1118 a_i(ji,jj,klbnd) = v_i(ji,jj,klbnd) / hi_max_typ(0,ntyp) 1119 ht_i(ji,jj,klbnd) = hi_max_typ(0,ntyp) 1103 DO ji = 1, jpi 1104 1105 IF (a_i(ji,jj,klbnd) > zeps) THEN 1106 IF (ht_i(ji,jj,klbnd) .LE. hi_max_typ(0,ntyp) .AND. hi_max_typ(0,ntyp) .GT. 0.0 ) THEN 1107 a_i(ji,jj,klbnd) = v_i(ji,jj,klbnd) / hi_max_typ(0,ntyp) 1108 ht_i(ji,jj,klbnd) = hi_max_typ(0,ntyp) 1109 ENDIF 1120 1110 ENDIF 1121 ENDIF 1122 END DO ! i 1111 END DO ! i 1123 1112 END DO ! j 1124 1113 1125 !------------------------------------------------------------------------------1126 ! 3) If a category thickness is not in bounds, shift the1127 ! entire area, volume, and energy to the neighboring category1128 !------------------------------------------------------------------------------1114 !------------------------------------------------------------------------------ 1115 ! 3) If a category thickness is not in bounds, shift the 1116 ! entire area, volume, and energy to the neighboring category 1117 !------------------------------------------------------------------------------ 1129 1118 !------------------------- 1130 1119 ! Initialize shift arrays … … 1133 1122 DO jl = klbnd, kubnd 1134 1123 DO jj = 1, jpj 1135 DO ji = 1, jpi1136 zdonor(ji,jj,jl) = 01137 zdaice(ji,jj,jl) = 0.01138 zdvice(ji,jj,jl) = 0.01139 END DO1124 DO ji = 1, jpi 1125 zdonor(ji,jj,jl) = 0 1126 zdaice(ji,jj,jl) = 0.0 1127 zdvice(ji,jj,jl) = 0.0 1128 END DO 1140 1129 END DO 1141 1130 END DO … … 1147 1136 DO jl = klbnd, kubnd - 1 ! loop over category boundaries 1148 1137 1149 !---------------------------------------1150 ! identify thicknesses that are too big1151 !---------------------------------------1138 !--------------------------------------- 1139 ! identify thicknesses that are too big 1140 !--------------------------------------- 1152 1141 zshiftflag = 0 1153 1142 … … 1166 1155 IF ( zshiftflag == 1 ) THEN 1167 1156 1168 !------------------------------1169 ! Shift ice between categories1170 !------------------------------1157 !------------------------------ 1158 ! Shift ice between categories 1159 !------------------------------ 1171 1160 CALL lim_itd_shiftice (klbnd, kubnd, zdonor, zdaice, zdvice) 1172 1173 !------------------------1174 ! Reset shift parameters1175 !------------------------1161 1162 !------------------------ 1163 ! Reset shift parameters 1164 !------------------------ 1176 1165 DO jj = 1, jpj 1177 DO ji = 1, jpi1178 zdonor(ji,jj,jl) = 01179 zdaice(ji,jj,jl) = 0.01180 zdvice(ji,jj,jl) = 0.01181 END DO1166 DO ji = 1, jpi 1167 zdonor(ji,jj,jl) = 0 1168 zdaice(ji,jj,jl) = 0.0 1169 zdvice(ji,jj,jl) = 0.0 1170 END DO 1182 1171 END DO 1183 1172 … … 1192 1181 DO jl = kubnd - 1, 1, -1 ! loop over category boundaries 1193 1182 1194 !-----------------------------------------1195 ! Identify thicknesses that are too small1196 !-----------------------------------------1183 !----------------------------------------- 1184 ! Identify thicknesses that are too small 1185 !----------------------------------------- 1197 1186 zshiftflag = 0 1198 1187 … … 1213 1202 IF (zshiftflag==1) THEN 1214 1203 1215 !------------------------------1216 ! Shift ice between categories1217 !------------------------------1204 !------------------------------ 1205 ! Shift ice between categories 1206 !------------------------------ 1218 1207 CALL lim_itd_shiftice (klbnd, kubnd, zdonor, zdaice, zdvice) 1219 1208 1220 !------------------------1221 ! Reset shift parameters1222 !------------------------1209 !------------------------ 1210 ! Reset shift parameters 1211 !------------------------ 1223 1212 DO jj = 1, jpj 1224 DO ji = 1, jpi 1225 zdonor(ji,jj,jl) = 0 1226 zdaice(ji,jj,jl) = 0.0 1227 zdvice(ji,jj,jl) = 0.0 1213 DO ji = 1, jpi 1214 zdonor(ji,jj,jl) = 0 1215 zdaice(ji,jj,jl) = 0.0 1216 zdvice(ji,jj,jl) = 0.0 1217 END DO 1228 1218 END DO 1229 END DO1230 1219 1231 1220 ENDIF ! zshiftflag … … 1233 1222 END DO ! jl 1234 1223 1235 !------------------------------------------------------------------------------1236 ! 4) Conservation check1237 !------------------------------------------------------------------------------1238 1239 IF ( con_i ) THEN1240 CALL lim_column_sum (jpl, v_i, vt_i_final)1241 fieldid = ' v_i : limitd_reb '1242 CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid)1243 1244 CALL lim_column_sum (jpl, v_s, vt_s_final)1245 fieldid = ' v_s : limitd_reb '1246 CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid)1247 ENDIF1248 1249 1224 !------------------------------------------------------------------------------ 1225 ! 4) Conservation check 1226 !------------------------------------------------------------------------------ 1227 1228 IF ( con_i ) THEN 1229 CALL lim_column_sum (jpl, v_i, vt_i_final) 1230 fieldid = ' v_i : limitd_reb ' 1231 CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) 1232 1233 CALL lim_column_sum (jpl, v_s, vt_s_final) 1234 fieldid = ' v_s : limitd_reb ' 1235 CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid) 1236 ENDIF 1237 1238 END SUBROUTINE lim_itd_th_reb 1250 1239 1251 1240 #else … … 1268 1257 END SUBROUTINE lim_itd_th_reb 1269 1258 #endif 1270 1259 END MODULE limitd_th
Note: See TracChangeset
for help on using the changeset viewer.
