Changeset 5965 for branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
- Timestamp:
- 2015-12-01T16:35:30+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
r4624 r5965 5 5 !!====================================================================== 6 6 !! History : LIM ! 2006-02 (M. Vancoppenolle) Original code 7 !! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in smsw & sfx_ mec7 !! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in smsw & sfx_dyn 8 8 !! 4.0 ! 2011-02 (G. Madec) dynamical allocation 9 9 !!---------------------------------------------------------------------- … … 18 18 USE thd_ice ! LIM thermodynamics 19 19 USE ice ! LIM variables 20 USE par_ice ! LIM parameters21 20 USE dom_ice ! LIM domain 22 USE limthd_lac ! LIM23 21 USE limvar ! LIM 24 USE limcons ! LIM25 USE in_out_manager ! I/O manager26 22 USE lbclnk ! lateral boundary condition - MPP exchanges 27 23 USE lib_mpp ! MPP library 28 24 USE wrk_nemo ! work arrays 29 25 USE prtctl ! Print control 30 ! Check budget (Rousset) 26 27 USE in_out_manager ! I/O manager 31 28 USE iom ! I/O manager 32 USE lib_fortran ! glob_sum29 USE lib_fortran ! glob_sum 33 30 USE limdiahsb 34 USE timing ! Timing 31 USE timing ! Timing 32 USE limcons ! conservation tests 35 33 36 34 IMPLICIT NONE … … 40 38 PUBLIC lim_itd_me_icestrength 41 39 PUBLIC lim_itd_me_init 42 PUBLIC lim_itd_me_zapsmall 43 PUBLIC lim_itd_me_alloc ! called by iceini.F90 44 45 REAL(wp) :: epsi20 = 1.e-20_wp ! constant values 46 REAL(wp) :: epsi10 = 1.e-10_wp ! constant values 47 REAL(wp) :: epsi06 = 1.e-06_wp ! constant values 40 PUBLIC lim_itd_me_alloc ! called by sbc_lim_init 48 41 49 42 !----------------------------------------------------------------------- … … 129 122 !! and Elizabeth C. Hunke, LANL are gratefully acknowledged 130 123 !!--------------------------------------------------------------------! 131 INTEGER :: ji, jj, jk, jl ! dummy loop index 132 INTEGER :: niter, nitermax = 20 ! local integer 133 LOGICAL :: asum_error ! flag for asum .ne. 1 124 INTEGER :: ji, jj, jk, jl ! dummy loop index 125 INTEGER :: niter ! local integer 134 126 INTEGER :: iterate_ridging ! if true, repeat the ridging 135 REAL(wp) :: w1, tmpfac! local scalar127 REAL(wp) :: za, zfac ! local scalar 136 128 CHARACTER (len = 15) :: fieldid 137 REAL(wp), POINTER, DIMENSION(:,:) :: closing_net ! net rate at which area is removed (1/s)138 ! (ridging ice area - area of new ridges) / dt139 REAL(wp), POINTER, DIMENSION(:,:) :: divu_adv ! divu as implied by transport scheme (1/s)140 REAL(wp), POINTER, DIMENSION(:,:) :: opning ! rate of opening due to divergence/shear141 REAL(wp), POINTER, DIMENSION(:,:) :: closing_gross ! rate at which area removed, not counting area of new ridges142 REAL(wp), POINTER, DIMENSION(:,:) :: msnow_mlt ! mass of snow added to ocean (kg m-2)143 REAL(wp), POINTER, DIMENSION(:,:) :: esnow_mlt ! energy needed to melt snow in ocean (J m-2)144 REAL(wp), POINTER, DIMENSION(:,:) :: vt_i_init, vt_i_final ! ice volume summed over categories145 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)146 REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset)147 ! mass and salt flux (clem)148 REAL(wp) , POINTER, DIMENSION(:,:,:) :: zviold, zvsold, zsmvold ! old ice volume...129 REAL(wp), POINTER, DIMENSION(:,:) :: closing_net ! net rate at which area is removed (1/s) 130 ! (ridging ice area - area of new ridges) / dt 131 REAL(wp), POINTER, DIMENSION(:,:) :: divu_adv ! divu as implied by transport scheme (1/s) 132 REAL(wp), POINTER, DIMENSION(:,:) :: opning ! rate of opening due to divergence/shear 133 REAL(wp), POINTER, DIMENSION(:,:) :: closing_gross ! rate at which area removed, not counting area of new ridges 134 REAL(wp), POINTER, DIMENSION(:,:) :: msnow_mlt ! mass of snow added to ocean (kg m-2) 135 REAL(wp), POINTER, DIMENSION(:,:) :: esnow_mlt ! energy needed to melt snow in ocean (J m-2) 136 REAL(wp), POINTER, DIMENSION(:,:) :: vt_i_init, vt_i_final ! ice volume summed over categories 137 ! 138 INTEGER, PARAMETER :: nitermax = 20 139 ! 140 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 149 141 !!----------------------------------------------------------------------------- 150 142 IF( nn_timing == 1 ) CALL timing_start('limitd_me') 151 143 152 CALL wrk_alloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final ) 153 154 CALL wrk_alloc( jpi, jpj, jpl, zviold, zvsold, zsmvold ) ! clem 155 156 IF( numit == nstart ) CALL lim_itd_me_init ! Initialization (first time-step only) 144 CALL wrk_alloc( jpi,jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final ) 157 145 158 146 IF(ln_ctl) THEN … … 162 150 163 151 IF( ln_limdyn ) THEN ! Start ridging and rafting ! 164 ! ------------------------------- 165 !- check conservation (C Rousset) 166 IF (ln_limdiahsb) THEN 167 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 168 zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 169 zchk_fw_b = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 170 zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 171 ENDIF 172 !- check conservation (C Rousset) 173 ! ------------------------------- 174 175 ! mass and salt flux init (clem) 176 zviold(:,:,:) = v_i(:,:,:) 177 zvsold(:,:,:) = v_s(:,:,:) 178 zsmvold(:,:,:) = smv_i(:,:,:) 152 153 ! conservation test 154 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 155 156 CALL lim_var_zapsmall 157 CALL lim_var_glo2eqv ! equivalent variables, requested for rafting 179 158 180 159 !-----------------------------------------------------------------------------! 181 160 ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons 182 161 !-----------------------------------------------------------------------------! 183 Cp = 0.5 * grav * (rau0-rhoic) * rhoic / rau0! proport const for PE162 Cp = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0 ! proport const for PE 184 163 ! 185 164 CALL lim_itd_me_ridgeprep ! prepare ridging … … 215 194 ! (thick, newly ridged ice). 216 195 217 closing_net(ji,jj) = Cs * 0.5 * ( Delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp )196 closing_net(ji,jj) = rn_cs * 0.5 * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp ) 218 197 219 198 ! 2.2 divu_adv … … 259 238 ! Reduce the closing rate if more than 100% of the open water 260 239 ! would be removed. Reduce the opening rate proportionately. 261 IF ( ato_i(ji,jj) .GT. epsi10 .AND. athorn(ji,jj,0) .GT. 0.0 ) THEN 262 w1 = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice 263 IF ( w1 .GT. ato_i(ji,jj)) THEN 264 tmpfac = ato_i(ji,jj) / w1 265 closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac 266 opning(ji,jj) = opning(ji,jj) * tmpfac 267 ENDIF !w1 268 ENDIF !at0i and athorn 269 270 END DO ! ji 271 END DO ! jj 240 za = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice 241 IF( za > epsi20 ) THEN 242 zfac = MIN( 1._wp, ato_i(ji,jj) / za ) 243 closing_gross(ji,jj) = closing_gross(ji,jj) * zfac 244 opning (ji,jj) = opning (ji,jj) * zfac 245 ENDIF 246 247 END DO 248 END DO 272 249 273 250 ! correction to closing rate / opening if excessive ice removal … … 275 252 ! Reduce the closing rate if more than 100% of any ice category 276 253 ! would be removed. Reduce the opening rate proportionately. 277 278 254 DO jl = 1, jpl 279 255 DO jj = 1, jpj 280 256 DO ji = 1, jpi 281 IF ( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp )THEN 282 w1 = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice 283 IF ( w1 > a_i(ji,jj,jl) ) THEN 284 tmpfac = a_i(ji,jj,jl) / w1 285 closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac 286 opning (ji,jj) = opning (ji,jj) * tmpfac 287 ENDIF 257 za = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice 258 IF( za > epsi20 ) THEN 259 zfac = MIN( 1._wp, a_i(ji,jj,jl) / za ) 260 closing_gross(ji,jj) = closing_gross(ji,jj) * zfac 261 opning (ji,jj) = opning (ji,jj) * zfac 288 262 ENDIF 289 END DO !ji290 END DO ! jj291 END DO !jl263 END DO 264 END DO 265 END DO 292 266 293 267 ! 3.3 Redistribute area, volume, and energy. … … 298 272 ! 3.4 Compute total area of ice plus open water after ridging. 299 273 !-----------------------------------------------------------------------------! 300 301 CALL lim_itd_me_asumr 274 ! This is in general not equal to one because of divergence during transport 275 asum(:,:) = ato_i(:,:) 276 DO jl = 1, jpl 277 asum(:,:) = asum(:,:) + a_i(:,:,jl) 278 END DO 302 279 303 280 ! 3.5 Do we keep on iterating ??? … … 310 287 DO jj = 1, jpj 311 288 DO ji = 1, jpi 312 IF (ABS(asum(ji,jj) - kamax ) .LT.epsi10) THEN289 IF (ABS(asum(ji,jj) - kamax ) < epsi10) THEN 313 290 closing_net(ji,jj) = 0._wp 314 291 opning (ji,jj) = 0._wp … … 346 323 ! Convert ridging rate diagnostics to correct units. 347 324 ! Update fresh water and heat fluxes due to snow melt. 348 349 asum_error = .false.350 351 325 DO jj = 1, jpj 352 326 DO ji = 1, jpi 353 354 IF(ABS(asum(ji,jj) - kamax) > epsi10 ) asum_error = .true.355 327 356 328 dardg1dt(ji,jj) = dardg1dt(ji,jj) * r1_rdtice … … 362 334 ! 5) Heat, salt and freshwater fluxes 363 335 !-----------------------------------------------------------------------------! 364 fmmec(ji,jj) = fmmec(ji,jj) + msnow_mlt(ji,jj) * r1_rdtice ! fresh water source for ocean365 fhmec(ji,jj) = fhmec(ji,jj) + esnow_mlt(ji,jj) * r1_rdtice ! heat sink for ocean336 wfx_snw(ji,jj) = wfx_snw(ji,jj) + msnow_mlt(ji,jj) * r1_rdtice ! fresh water source for ocean 337 hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + esnow_mlt(ji,jj) * r1_rdtice ! heat sink for ocean (<0, W.m-2) 366 338 367 339 END DO … … 369 341 370 342 ! Check if there is a ridging error 371 DO jj = 1, jpj 372 DO ji = 1, jpi 373 IF( ABS( asum(ji,jj) - kamax) > epsi10 ) THEN ! there is a bug 374 WRITE(numout,*) ' ' 375 WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj) 376 WRITE(numout,*) ' limitd_me ' 377 WRITE(numout,*) ' POINT : ', ji, jj 378 WRITE(numout,*) ' jpl, a_i, athorn ' 379 WRITE(numout,*) 0, ato_i(ji,jj), athorn(ji,jj,0) 380 DO jl = 1, jpl 381 WRITE(numout,*) jl, a_i(ji,jj,jl), athorn(ji,jj,jl) 382 END DO 383 ENDIF ! asum 384 385 END DO !ji 386 END DO !jj 343 IF( lwp ) THEN 344 DO jj = 1, jpj 345 DO ji = 1, jpi 346 IF( ABS( asum(ji,jj) - kamax) > epsi10 ) THEN ! there is a bug 347 WRITE(numout,*) ' ' 348 WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj) 349 WRITE(numout,*) ' limitd_me ' 350 WRITE(numout,*) ' POINT : ', ji, jj 351 WRITE(numout,*) ' jpl, a_i, athorn ' 352 WRITE(numout,*) 0, ato_i(ji,jj), athorn(ji,jj,0) 353 DO jl = 1, jpl 354 WRITE(numout,*) jl, a_i(ji,jj,jl), athorn(ji,jj,jl) 355 END DO 356 ENDIF 357 END DO 358 END DO 359 END IF 387 360 388 361 ! Conservation check … … 393 366 ENDIF 394 367 368 CALL lim_var_agg( 1 ) 369 395 370 !-----------------------------------------------------------------------------! 396 ! 6) Updating state variables and trend terms (done in limupdate)371 ! control prints 397 372 !-----------------------------------------------------------------------------! 398 CALL lim_var_glo2eqv 399 CALL lim_itd_me_zapsmall 400 401 !-------------------------------- 402 ! Update mass/salt fluxes (clem) 403 !-------------------------------- 404 DO jl = 1, jpl 405 DO jj = 1, jpj 406 DO ji = 1, jpi 407 diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * r1_rdtice 408 rdm_ice(ji,jj) = rdm_ice(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * rhoic 409 rdm_snw(ji,jj) = rdm_snw(ji,jj) + ( v_s(ji,jj,jl) - zvsold(ji,jj,jl) ) * rhosn 410 sfx_mec(ji,jj) = sfx_mec(ji,jj) - ( smv_i(ji,jj,jl) - zsmvold(ji,jj,jl) ) * rhoic * r1_rdtice 411 END DO 412 END DO 413 END DO 414 415 IF(ln_ctl) THEN ! Control print 373 IF(ln_ctl) THEN 374 CALL lim_var_glo2eqv 375 416 376 CALL prt_ctl_info(' ') 417 377 CALL prt_ctl_info(' - Cell values : ') 418 378 CALL prt_ctl_info(' ~~~~~~~~~~~~~ ') 419 CALL prt_ctl(tab2d_1= area, clinfo1=' lim_itd_me : cell area :')379 CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_itd_me : cell area :') 420 380 CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_me : at_i :') 421 381 CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_me : vt_i :') … … 445 405 ENDIF 446 406 447 ! ------------------------------- 448 !- check conservation (C Rousset) 449 IF (ln_limdiahsb) THEN 450 zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 451 zchk_fw = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 452 453 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice 454 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 455 456 zchk_vmin = glob_min(v_i) 457 zchk_amax = glob_max(SUM(a_i,dim=3)) 458 zchk_amin = glob_min(a_i) 459 460 IF(lwp) THEN 461 IF ( ABS( zchk_v_i ) > 1.e-5 ) WRITE(numout,*) 'violation volume [m3/day] (limitd_me) = ',(zchk_v_i * rday) 462 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_me) = ',(zchk_smv * rday) 463 IF ( zchk_vmin < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limitd_me) = ',(zchk_vmin * 1.e-3) 464 IF ( zchk_amax > kamax+epsi10 ) WRITE(numout,*) 'violation a_i>amax (limitd_me) = ',zchk_amax 465 IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limitd_me) = ',zchk_amin 466 ENDIF 467 ENDIF 468 !- check conservation (C Rousset) 469 ! ------------------------------- 407 ! conservation test 408 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 470 409 471 410 ENDIF ! ln_limdyn=.true. 472 411 ! 473 412 CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final ) 474 !475 CALL wrk_dealloc( jpi, jpj, jpl, zviold, zvsold, zsmvold ) ! clem476 413 ! 477 414 IF( nn_timing == 1 ) CALL timing_stop('limitd_me') … … 494 431 !!---------------------------------------------------------------------- 495 432 INTEGER, INTENT(in) :: kstrngth ! = 1 for Rothrock formulation, 0 for Hibler (1979) 496 497 INTEGER :: ji,jj, jl ! dummy loop indices 498 INTEGER :: ksmooth ! smoothing the resistance to deformation 499 INTEGER :: numts_rm ! number of time steps for the P smoothing 500 REAL(wp) :: hi, zw1, zp, zdummy, zzc, z1_3 ! local scalars 433 INTEGER :: ji,jj, jl ! dummy loop indices 434 INTEGER :: ksmooth ! smoothing the resistance to deformation 435 INTEGER :: numts_rm ! number of time steps for the P smoothing 436 REAL(wp) :: zhi, zp, z1_3 ! local scalars 501 437 REAL(wp), POINTER, DIMENSION(:,:) :: zworka ! temporary array used here 502 438 !!---------------------------------------------------------------------- … … 524 460 ! 525 461 IF( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp ) THEN 526 hi = v_i(ji,jj,jl) / a_i(ji,jj,jl)462 zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 527 463 !---------------------------- 528 464 ! PE loss from deforming ice 529 465 !---------------------------- 530 strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * hi *hi466 strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * zhi * zhi 531 467 532 468 !-------------------------- 533 469 ! PE gain from rafting ice 534 470 !-------------------------- 535 strength(ji,jj) = strength(ji,jj) + 2._wp * araft(ji,jj,jl) * hi *hi471 strength(ji,jj) = strength(ji,jj) + 2._wp * araft(ji,jj,jl) * zhi * zhi 536 472 537 473 !---------------------------- 538 474 ! PE gain from ridging ice 539 475 !---------------------------- 540 strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl) /krdg(ji,jj,jl) &541 * z1_3 * ( hrmax(ji,jj,jl)**3 - hrmin(ji,jj,jl)**3) / ( hrmax(ji,jj,jl)-hrmin(ji,jj,jl) )542 !!gm Optimization: (a**3-b**3)/(a-b) = a*a+ab+b*b ==> less costly operations even if a**3 is replaced by a*a*a...543 ENDIF ! aicen > epsi10476 strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl) / krdg(ji,jj,jl) & 477 * z1_3 * ( hrmax(ji,jj,jl)**2 + hrmin(ji,jj,jl)**2 + hrmax(ji,jj,jl) * hrmin(ji,jj,jl) ) 478 !!(a**3-b**3)/(a-b) = a*a+ab+b*b 479 ENDIF 544 480 ! 545 END DO ! ji 546 END DO !jj 547 END DO !jl 548 549 zzc = Cf * Cp ! where Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) and Cf accounts for frictional dissipation 550 strength(:,:) = zzc * strength(:,:) / aksum(:,:) 551 481 END DO 482 END DO 483 END DO 484 485 strength(:,:) = rn_pe_rdg * Cp * strength(:,:) / aksum(:,:) 486 ! where Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation 552 487 ksmooth = 1 553 488 … … 557 492 ELSE ! kstrngth ne 1: Hibler (1979) form 558 493 ! 559 strength(:,:) = Pstar * vt_i(:,:) * EXP( - C_rhg * ( 1._wp - at_i(:,:) ) )494 strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) ) ) 560 495 ! 561 496 ksmooth = 1 … … 569 504 ! CAN BE REMOVED 570 505 ! 571 IF ( brinstren_swi == 1) THEN506 IF( ln_icestr_bvf ) THEN 572 507 573 508 DO jj = 1, jpj 574 509 DO ji = 1, jpi 575 IF ( bv_i(ji,jj) .GT. 0.0 ) THEN576 zdummy = MIN ( bv_i(ji,jj), 0.10 ) * MIN( bv_i(ji,jj), 0.10 )577 ELSE578 zdummy = 0.0579 ENDIF580 510 strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bv_i(ji,jj),0.0))) 581 END DO ! j582 END DO ! i511 END DO 512 END DO 583 513 584 514 ENDIF … … 596 526 CALL lbc_lnk( strength, 'T', 1. ) 597 527 598 DO jj = 2, jpj - 1 599 DO ji = 2, jpi - 1 600 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi10) THEN ! ice is 601 ! present 602 zworka(ji,jj) = 4.0 * strength(ji,jj) & 603 & + strength(ji-1,jj) * tms(ji-1,jj) & 604 & + strength(ji+1,jj) * tms(ji+1,jj) & 605 & + strength(ji,jj-1) * tms(ji,jj-1) & 606 & + strength(ji,jj+1) * tms(ji,jj+1) 607 608 zw1 = 4.0 + tms(ji-1,jj) + tms(ji+1,jj) + tms(ji,jj-1) + tms(ji,jj+1) 609 zworka(ji,jj) = zworka(ji,jj) / zw1 528 DO jj = 2, jpjm1 529 DO ji = 2, jpim1 530 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp) THEN 531 zworka(ji,jj) = ( 4.0 * strength(ji,jj) & 532 & + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) & 533 & + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) & 534 & ) / ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) ) 610 535 ELSE 611 536 zworka(ji,jj) = 0._wp … … 614 539 END DO 615 540 616 DO jj = 2, jpj -1617 DO ji = 2, jpi -1541 DO jj = 2, jpjm1 542 DO ji = 2, jpim1 618 543 strength(ji,jj) = zworka(ji,jj) 619 544 END DO … … 621 546 CALL lbc_lnk( strength, 'T', 1. ) 622 547 623 ENDIF ! ksmooth548 ENDIF 624 549 625 550 !-------------------- … … 638 563 DO jj = 1, jpj - 1 639 564 DO ji = 1, jpi - 1 640 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi10) THEN ! ice is present565 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp) THEN 641 566 numts_rm = 1 ! number of time steps for the running mean 642 IF ( strp1(ji,jj) .GT.0.0 ) numts_rm = numts_rm + 1643 IF ( strp2(ji,jj) .GT.0.0 ) numts_rm = numts_rm + 1567 IF ( strp1(ji,jj) > 0.0 ) numts_rm = numts_rm + 1 568 IF ( strp2(ji,jj) > 0.0 ) numts_rm = numts_rm + 1 644 569 zp = ( strength(ji,jj) + strp1(ji,jj) + strp2(ji,jj) ) / numts_rm 645 570 strp2(ji,jj) = strp1(ji,jj) … … 670 595 !!---------------------------------------------------------------------! 671 596 INTEGER :: ji,jj, jl ! dummy loop indices 672 INTEGER :: krdg_index ! 673 REAL(wp) :: Gstari, astari, hi, hrmean, zdummy ! local scalar 597 REAL(wp) :: Gstari, astari, zhi, hrmean, zdummy ! local scalar 674 598 REAL(wp), POINTER, DIMENSION(:,:) :: zworka ! temporary array used here 675 599 REAL(wp), POINTER, DIMENSION(:,:,:) :: Gsum ! Gsum(n) = sum of areas in categories 0 to n … … 679 603 CALL wrk_alloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 ) 680 604 681 Gstari = 1.0/ Gstar682 astari = 1.0/ astar605 Gstari = 1.0/rn_gstar 606 astari = 1.0/rn_astar 683 607 aksum(:,:) = 0.0 684 608 athorn(:,:,:) = 0.0 … … 691 615 692 616 ! ! Zero out categories with very small areas 693 CALL lim_ itd_me_zapsmall617 CALL lim_var_zapsmall 694 618 695 619 !------------------------------------------------------------------------------! … … 698 622 699 623 ! Compute total area of ice plus open water. 700 CALL lim_itd_me_asumr 701 ! This is in general not equal to one 702 ! because of divergence during transport 624 ! This is in general not equal to one because of divergence during transport 625 asum(:,:) = ato_i(:,:) 626 DO jl = 1, jpl 627 asum(:,:) = asum(:,:) + a_i(:,:,jl) 628 END DO 703 629 704 630 ! Compute cumulative thickness distribution function … … 708 634 709 635 Gsum(:,:,-1) = 0._wp 710 711 DO jj = 1, jpj 712 DO ji = 1, jpi 713 IF( ato_i(ji,jj) > epsi10 ) THEN ; Gsum(ji,jj,0) = ato_i(ji,jj) 714 ELSE ; Gsum(ji,jj,0) = 0._wp 715 ENDIF 716 END DO 717 END DO 636 Gsum(:,:,0 ) = ato_i(:,:) 718 637 719 638 ! for each value of h, you have to add ice concentration then 720 639 DO jl = 1, jpl 721 DO jj = 1, jpj 722 DO ji = 1, jpi 723 IF( a_i(ji,jj,jl) .GT. epsi10 ) THEN ; Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl) 724 ELSE ; Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) 725 ENDIF 726 END DO 727 END DO 640 Gsum(:,:,jl) = Gsum(:,:,jl-1) + a_i(:,:,jl) 728 641 END DO 729 642 … … 746 659 !----------------------------------------------------------------- 747 660 748 krdg_index = 1 749 750 IF( krdg_index == 0 ) THEN !--- Linear formulation (Thorndike et al., 1975) 751 DO jl = 0, ice_cat_bounds(1,2) ! only undeformed ice participates 661 IF( nn_partfun == 0 ) THEN !--- Linear formulation (Thorndike et al., 1975) 662 DO jl = 0, jpl 752 663 DO jj = 1, jpj 753 664 DO ji = 1, jpi 754 IF( Gsum(ji,jj,jl) < Gstar) THEN755 athorn(ji,jj,jl) = Gstari * ( Gsum(ji,jj,jl)-Gsum(ji,jj,jl-1)) * &756 (2.0 - (Gsum(ji,jj,jl-1)+Gsum(ji,jj,jl))*Gstari)757 ELSEIF (Gsum(ji,jj,jl-1) < Gstar) THEN758 athorn(ji,jj,jl) = Gstari * ( Gstar-Gsum(ji,jj,jl-1)) * &759 (2.0 - (Gsum(ji,jj,jl-1)+Gstar)*Gstari)665 IF( Gsum(ji,jj,jl) < rn_gstar) THEN 666 athorn(ji,jj,jl) = Gstari * ( Gsum(ji,jj,jl) - Gsum(ji,jj,jl-1) ) * & 667 & ( 2.0 - (Gsum(ji,jj,jl-1) + Gsum(ji,jj,jl) ) * Gstari ) 668 ELSEIF (Gsum(ji,jj,jl-1) < rn_gstar) THEN 669 athorn(ji,jj,jl) = Gstari * ( rn_gstar - Gsum(ji,jj,jl-1) ) * & 670 & ( 2.0 - ( Gsum(ji,jj,jl-1) + rn_gstar ) * Gstari ) 760 671 ELSE 761 672 athorn(ji,jj,jl) = 0.0 762 673 ENDIF 763 END DO ! ji764 END DO ! jj765 END DO ! jl674 END DO 675 END DO 676 END DO 766 677 767 678 ELSE !--- Exponential, more stable formulation (Lipscomb et al, 2007) 768 679 ! 769 680 zdummy = 1._wp / ( 1._wp - EXP(-astari) ) ! precompute exponential terms using Gsum as a work array 770 771 681 DO jl = -1, jpl 772 682 Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy 773 END DO !jl774 DO jl = 0, ice_cat_bounds(1,2)683 END DO 684 DO jl = 0, jpl 775 685 athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl) 776 686 END DO 777 687 ! 778 ENDIF ! krdg_index779 780 IF( raftswi == 1) THEN ! Ridging and rafting ice participation functions688 ENDIF 689 690 IF( ln_rafting ) THEN ! Ridging and rafting ice participation functions 781 691 ! 782 692 DO jl = 1, jpl 783 693 DO jj = 1, jpj 784 694 DO ji = 1, jpi 785 IF ( athorn(ji,jj,jl) .GT.0._wp ) THEN695 IF ( athorn(ji,jj,jl) > 0._wp ) THEN 786 696 !!gm TANH( -X ) = - TANH( X ) so can be computed only 1 time.... 787 aridge(ji,jj,jl) = ( TANH ( Craft * ( ht_i(ji,jj,jl) - hparmeter) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl)788 araft (ji,jj,jl) = ( TANH ( - Craft * ( ht_i(ji,jj,jl) - hparmeter) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl)697 aridge(ji,jj,jl) = ( TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl) 698 araft (ji,jj,jl) = ( TANH ( -rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl) 789 699 IF ( araft(ji,jj,jl) < epsi06 ) araft(ji,jj,jl) = 0._wp 790 700 aridge(ji,jj,jl) = MAX( athorn(ji,jj,jl) - araft(ji,jj,jl), 0.0 ) 791 ENDIF ! athorn792 END DO ! ji793 END DO ! jj794 END DO ! jl795 796 ELSE ! raftswi = 0701 ENDIF 702 END DO 703 END DO 704 END DO 705 706 ELSE 797 707 ! 798 708 DO jl = 1, jpl … … 802 712 ENDIF 803 713 804 IF ( raftswi == 1) THEN805 806 IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) .GT. epsi10) THEN714 IF( ln_rafting ) THEN 715 716 IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) > epsi10 .AND. lwp ) THEN 807 717 DO jl = 1, jpl 808 718 DO jj = 1, jpj 809 719 DO ji = 1, jpi 810 IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) .GT.epsi10 ) THEN720 IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) > epsi10 ) THEN 811 721 WRITE(numout,*) ' ALERTE 96 : wrong participation function ... ' 812 722 WRITE(numout,*) ' ji, jj, jl : ', ji, jj, jl … … 854 764 DO ji = 1, jpi 855 765 856 IF (a_i(ji,jj,jl) .GT. epsi10 .AND. athorn(ji,jj,jl) .GT.0.0 ) THEN857 hi = v_i(ji,jj,jl) / a_i(ji,jj,jl)858 hrmean = MAX(SQRT( Hstar*hi),hi*krdgmin)859 hrmin(ji,jj,jl) = MIN(2.0* hi, 0.5*(hrmean +hi))766 IF (a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0.0 ) THEN 767 zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 768 hrmean = MAX(SQRT(rn_hstar*zhi), zhi*krdgmin) 769 hrmin(ji,jj,jl) = MIN(2.0*zhi, 0.5*(hrmean + zhi)) 860 770 hrmax(ji,jj,jl) = 2.0*hrmean - hrmin(ji,jj,jl) 861 hraft(ji,jj,jl) = kraft* hi862 krdg(ji,jj,jl) = hrmean / hi771 hraft(ji,jj,jl) = kraft*zhi 772 krdg(ji,jj,jl) = hrmean / zhi 863 773 ELSE 864 774 hraft(ji,jj,jl) = 0.0 … … 868 778 ENDIF 869 779 870 END DO ! ji871 END DO ! jj872 END DO ! jl780 END DO 781 END DO 782 END DO 873 783 874 784 ! Normalization factor : aksum, ensures mass conservation … … 902 812 LOGICAL, PARAMETER :: l_conservation_check = .true. ! if true, check conservation (useful for debugging) 903 813 ! 904 LOGICAL :: neg_ato_i ! flag for ato_i(i,j) < -puny905 LOGICAL :: large_afrac ! flag for afrac > 1906 LOGICAL :: large_afrft ! flag for afrac > 1907 814 INTEGER :: ji, jj, jl, jl1, jl2, jk ! dummy loop indices 908 815 INTEGER :: ij ! horizontal index, combines i and j loops 909 816 INTEGER :: icells ! number of cells with aicen > puny 910 REAL(wp) :: zindb, zsrdg2 ! local scalar 911 REAL(wp) :: hL, hR, farea, zdummy, zdummy0, ztmelts ! left and right limits of integration 817 REAL(wp) :: hL, hR, farea, ztmelts ! left and right limits of integration 912 818 913 819 INTEGER , POINTER, DIMENSION(:) :: indxi, indxj ! compressed indices … … 917 823 918 824 REAL(wp), POINTER, DIMENSION(:,:,:) :: aicen_init, vicen_init ! ice area & volume before ridging 919 REAL(wp), POINTER, DIMENSION(:,:,:) :: vsn on_init, esnon_init ! snow volume & energy before ridging825 REAL(wp), POINTER, DIMENSION(:,:,:) :: vsnwn_init, esnwn_init ! snow volume & energy before ridging 920 826 REAL(wp), POINTER, DIMENSION(:,:,:) :: smv_i_init, oa_i_init ! ice salinity & age before ridging 921 827 … … 925 831 REAL(wp), POINTER, DIMENSION(:,:) :: ardg1 , ardg2 ! area of ice ridged & new ridges 926 832 REAL(wp), POINTER, DIMENSION(:,:) :: vsrdg , esrdg ! snow volume & energy of ridging ice 927 REAL(wp), POINTER, DIMENSION(:,:) :: oirdg1, oirdg2 ! areal age content of ridged & rifging ice928 833 REAL(wp), POINTER, DIMENSION(:,:) :: dhr , dhr2 ! hrmax - hrmin & hrmax^2 - hrmin^2 929 834 … … 934 839 REAL(wp), POINTER, DIMENSION(:,:) :: srdg2 ! sal*volume of new ridges 935 840 REAL(wp), POINTER, DIMENSION(:,:) :: smsw ! sal*volume of water trapped into ridges 841 REAL(wp), POINTER, DIMENSION(:,:) :: oirdg1, oirdg2 ! ice age of ice ridged 936 842 937 843 REAL(wp), POINTER, DIMENSION(:,:) :: afrft ! fraction of category area rafted … … 939 845 REAL(wp), POINTER, DIMENSION(:,:) :: virft , vsrft ! ice & snow volume of rafting ice 940 846 REAL(wp), POINTER, DIMENSION(:,:) :: esrft , smrft ! snow energy & salinity of rafting ice 941 REAL(wp), POINTER, DIMENSION(:,:) :: oirft1, oirft2 ! areal age content of rafted ice & rafting ice847 REAL(wp), POINTER, DIMENSION(:,:) :: oirft1, oirft2 ! ice age of ice rafted 942 848 943 849 REAL(wp), POINTER, DIMENSION(:,:,:) :: eirft ! ice energy of rafting ice … … 947 853 !!---------------------------------------------------------------------- 948 854 949 CALL wrk_alloc( (jpi+1)*(jpj+1), indxi, indxj )950 CALL wrk_alloc( jpi, jpj, vice_init, vice_final, eice_init, eice_final )951 CALL wrk_alloc( jpi, jpj, afrac, fvol , ardg1, ardg2, vsrdg, esrdg, oirdg1, oirdg2, dhr, dhr2 )952 CALL wrk_alloc( jpi, jpj, vrdg1, vrdg2, vsw , srdg1, srdg2, smsw)953 CALL wrk_alloc( jpi, jpj, afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )954 CALL wrk_alloc( jpi, jpj, jpl, aicen_init, vicen_init, vsnon_init, esnon_init, smv_i_init, oa_i_init )955 CALL wrk_alloc( jpi, jpj, jkmax, eirft, erdg1, erdg2, ersw )956 CALL wrk_alloc( jpi, jpj, jkmax, jpl, eicen_init )855 CALL wrk_alloc( (jpi+1)*(jpj+1), indxi, indxj ) 856 CALL wrk_alloc( jpi, jpj, vice_init, vice_final, eice_init, eice_final ) 857 CALL wrk_alloc( jpi, jpj, afrac, fvol , ardg1, ardg2, vsrdg, esrdg, dhr, dhr2 ) 858 CALL wrk_alloc( jpi, jpj, vrdg1, vrdg2, vsw , srdg1, srdg2, smsw, oirdg1, oirdg2 ) 859 CALL wrk_alloc( jpi, jpj, afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 ) 860 CALL wrk_alloc( jpi, jpj, jpl, aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init ) 861 CALL wrk_alloc( jpi, jpj, nlay_i, eirft, erdg1, erdg2, ersw ) 862 CALL wrk_alloc( jpi, jpj, nlay_i, jpl, eicen_init ) 957 863 958 864 ! Conservation check … … 962 868 CALL lim_column_sum (jpl, v_i, vice_init ) 963 869 CALL lim_column_sum_energy (jpl, nlay_i, e_i, eice_init ) 964 DO ji = mi0( jiindx), mi1(jiindx)965 DO jj = mj0(j jindx), mj1(jjindx)870 DO ji = mi0(iiceprt), mi1(iiceprt) 871 DO jj = mj0(jiceprt), mj1(jiceprt) 966 872 WRITE(numout,*) ' vice_init : ', vice_init(ji,jj) 967 873 WRITE(numout,*) ' eice_init : ', eice_init(ji,jj) … … 973 879 ! 1) Compute change in open water area due to closing and opening. 974 880 !------------------------------------------------------------------------------- 975 976 neg_ato_i = .false.977 978 881 DO jj = 1, jpj 979 882 DO ji = 1, jpi 980 883 ato_i(ji,jj) = ato_i(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice & 981 884 & + opning(ji,jj) * rdt_ice 982 IF ( ato_i(ji,jj) < -epsi10 ) THEN983 neg_ato_i = .TRUE.984 ELSEIF( ato_i(ji,jj) < 0._wp ) THEN ! roundoff error885 IF ( ato_i(ji,jj) < -epsi10 ) THEN ! there is a bug 886 IF(lwp) WRITE(numout,*) 'Ridging error: ato_i < 0 -- ato_i : ',ato_i(ji,jj) 887 ELSEIF( ato_i(ji,jj) < 0._wp ) THEN ! roundoff error 985 888 ato_i(ji,jj) = 0._wp 986 889 ENDIF 987 END DO !jj 988 END DO !ji 989 990 ! if negative open water area alert it 991 IF( neg_ato_i ) THEN ! there is a bug 992 DO jj = 1, jpj 993 DO ji = 1, jpi 994 IF( ato_i(ji,jj) < -epsi10 ) THEN 995 WRITE(numout,*) '' 996 WRITE(numout,*) 'Ridging error: ato_i < 0' 997 WRITE(numout,*) 'ato_i : ', ato_i(ji,jj) 998 ENDIF ! ato_i < -epsi10 999 END DO 1000 END DO 1001 ENDIF 890 END DO 891 END DO 1002 892 1003 893 !----------------------------------------------------------------- 1004 894 ! 2) Save initial state variables 1005 895 !----------------------------------------------------------------- 1006 1007 DO jl = 1, jpl 1008 aicen_init(:,:,jl) = a_i(:,:,jl) 1009 vicen_init(:,:,jl) = v_i(:,:,jl) 1010 vsnon_init(:,:,jl) = v_s(:,:,jl) 1011 ! 1012 smv_i_init(:,:,jl) = smv_i(:,:,jl) 1013 oa_i_init (:,:,jl) = oa_i (:,:,jl) 1014 END DO !jl 1015 1016 esnon_init(:,:,:) = e_s(:,:,1,:) 1017 1018 DO jl = 1, jpl 1019 DO jk = 1, nlay_i 1020 eicen_init(:,:,jk,jl) = e_i(:,:,jk,jl) 1021 END DO 1022 END DO 896 aicen_init(:,:,:) = a_i (:,:,:) 897 vicen_init(:,:,:) = v_i (:,:,:) 898 vsnwn_init(:,:,:) = v_s (:,:,:) 899 smv_i_init(:,:,:) = smv_i(:,:,:) 900 esnwn_init(:,:,:) = e_s (:,:,1,:) 901 eicen_init(:,:,:,:) = e_i (:,:,:,:) 902 oa_i_init (:,:,:) = oa_i (:,:,:) 1023 903 1024 904 ! … … 1043 923 indxi(icells) = ji 1044 924 indxj(icells) = jj 1045 ENDIF ! test on a_icen_init 1046 END DO ! ji 1047 END DO ! jj 1048 1049 large_afrac = .false. 1050 large_afrft = .false. 1051 1052 !CDIR NODEP 925 ENDIF 926 END DO 927 END DO 928 1053 929 DO ij = 1, icells 1054 930 ji = indxi(ij) … … 1064 940 arft2(ji,jj) = arft1(ji,jj) / kraft 1065 941 1066 oirdg1(ji,jj)= aridge(ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice1067 oirft1(ji,jj)= araft (ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice1068 oirdg2(ji,jj)= oirdg1(ji,jj) / krdg(ji,jj,jl1)1069 oirft2(ji,jj)= oirft1(ji,jj) / kraft1070 1071 942 !--------------------------------------------------------------- 1072 943 ! 3.3) Compute ridging /rafting fractions, make sure afrac <=1 … … 1076 947 afrft(ji,jj) = arft1(ji,jj) / aicen_init(ji,jj,jl1) !rafting 1077 948 1078 IF (afrac(ji,jj) > kamax + epsi10) THEN !riging1079 large_afrac = .true.1080 ELSEIF (afrac(ji,jj) > kamax) THEN! roundoff error949 IF( afrac(ji,jj) > kamax + epsi10 ) THEN ! there is a bug 950 IF(lwp) WRITE(numout,*) ' ardg > a_i -- ardg, aicen_init : ', ardg1(ji,jj), aicen_init(ji,jj,jl1) 951 ELSEIF( afrac(ji,jj) > kamax ) THEN ! roundoff error 1081 952 afrac(ji,jj) = kamax 1082 953 ENDIF 1083 IF (afrft(ji,jj) > kamax + epsi10) THEN !rafting 1084 large_afrft = .true. 1085 ELSEIF (afrft(ji,jj) > kamax) THEN ! roundoff error 954 955 IF( afrft(ji,jj) > kamax + epsi10 ) THEN ! there is a bug 956 IF(lwp) WRITE(numout,*) ' arft > a_i -- arft, aicen_init : ', arft1(ji,jj), aicen_init(ji,jj,jl1) 957 ELSEIF( afrft(ji,jj) > kamax) THEN ! roundoff error 1086 958 afrft(ji,jj) = kamax 1087 959 ENDIF … … 1091 963 ! / rafting category n1. 1092 964 !-------------------------------------------------------------------------- 1093 vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj) / ( 1._wp + ridge_por ) 1094 vrdg2(ji,jj) = vrdg1(ji,jj) * ( 1. + ridge_por ) 1095 vsw (ji,jj) = vrdg1(ji,jj) * ridge_por 1096 1097 vsrdg(ji,jj) = vsnon_init(ji,jj,jl1) * afrac(ji,jj) 1098 esrdg(ji,jj) = esnon_init(ji,jj,jl1) * afrac(ji,jj) 1099 srdg1(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) / ( 1._wp + ridge_por ) 1100 srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) 965 vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj) 966 vrdg2(ji,jj) = vrdg1(ji,jj) * ( 1. + rn_por_rdg ) 967 vsw (ji,jj) = vrdg1(ji,jj) * rn_por_rdg 968 969 vsrdg (ji,jj) = vsnwn_init(ji,jj,jl1) * afrac(ji,jj) 970 esrdg (ji,jj) = esnwn_init(ji,jj,jl1) * afrac(ji,jj) 971 srdg1 (ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) 972 oirdg1(ji,jj) = oa_i_init (ji,jj,jl1) * afrac(ji,jj) 973 oirdg2(ji,jj) = oa_i_init (ji,jj,jl1) * afrac(ji,jj) / krdg(ji,jj,jl1) 1101 974 1102 975 ! rafting volumes, heat contents ... 1103 virft(ji,jj) = vicen_init(ji,jj,jl1) * afrft(ji,jj) 1104 vsrft(ji,jj) = vsnon_init(ji,jj,jl1) * afrft(ji,jj) 1105 esrft(ji,jj) = esnon_init(ji,jj,jl1) * afrft(ji,jj) 1106 smrft(ji,jj) = smv_i_init(ji,jj,jl1) * afrft(ji,jj) 976 virft (ji,jj) = vicen_init(ji,jj,jl1) * afrft(ji,jj) 977 vsrft (ji,jj) = vsnwn_init(ji,jj,jl1) * afrft(ji,jj) 978 esrft (ji,jj) = esnwn_init(ji,jj,jl1) * afrft(ji,jj) 979 smrft (ji,jj) = smv_i_init(ji,jj,jl1) * afrft(ji,jj) 980 oirft1(ji,jj) = oa_i_init (ji,jj,jl1) * afrft(ji,jj) 981 oirft2(ji,jj) = oa_i_init (ji,jj,jl1) * afrft(ji,jj) / kraft 1107 982 1108 983 ! substract everything 1109 a_i(ji,jj,jl1) = a_i(ji,jj,jl1) - ardg1(ji,jj) - arft1(ji,jj) 1110 v_i(ji,jj,jl1) = v_i(ji,jj,jl1) - vrdg1(ji,jj) - virft(ji,jj) 1111 v_s(ji,jj,jl1) = v_s(ji,jj,jl1) - vsrdg(ji,jj) - vsrft(ji,jj) 1112 e_s(ji,jj,1,jl1) = e_s(ji,jj,1,jl1) - esrdg(ji,jj) - esrft(ji,jj) 984 a_i(ji,jj,jl1) = a_i(ji,jj,jl1) - ardg1 (ji,jj) - arft1 (ji,jj) 985 v_i(ji,jj,jl1) = v_i(ji,jj,jl1) - vrdg1 (ji,jj) - virft (ji,jj) 986 v_s(ji,jj,jl1) = v_s(ji,jj,jl1) - vsrdg (ji,jj) - vsrft (ji,jj) 987 e_s(ji,jj,1,jl1) = e_s(ji,jj,1,jl1) - esrdg (ji,jj) - esrft (ji,jj) 988 smv_i(ji,jj,jl1) = smv_i(ji,jj,jl1) - srdg1 (ji,jj) - smrft (ji,jj) 1113 989 oa_i(ji,jj,jl1) = oa_i(ji,jj,jl1) - oirdg1(ji,jj) - oirft1(ji,jj) 1114 smv_i(ji,jj,jl1) = smv_i(ji,jj,jl1) - srdg1(ji,jj) - smrft(ji,jj)1115 990 1116 991 !----------------------------------------------------------------- 1117 992 ! 3.5) Compute properties of new ridges 1118 993 !----------------------------------------------------------------- 1119 !--------- ----994 !--------- 1120 995 ! Salinity 1121 !------------- 1122 smsw(ji,jj) = sss_m(ji,jj) * vsw(ji,jj) * rhoic / rau0 ! salt content of seawater frozen in voids 1123 1124 zsrdg2 = srdg1(ji,jj) + smsw(ji,jj) ! salt content of new ridge 1125 1126 srdg2(ji,jj) = MIN( s_i_max * vrdg2(ji,jj) , zsrdg2 ) ! impose a maximum salinity 996 !--------- 997 smsw(ji,jj) = vsw(ji,jj) * sss_m(ji,jj) ! salt content of seawater frozen in voids !! MV HC2014 998 srdg2(ji,jj) = srdg1(ji,jj) + smsw(ji,jj) ! salt content of new ridge 999 1000 !srdg2(ji,jj) = MIN( rn_simax * vrdg2(ji,jj) , zsrdg2 ) ! impose a maximum salinity 1127 1001 1128 ! ! excess of salt is flushed into the ocean 1129 !sfx_mec(ji,jj) = sfx_mec(ji,jj) + ( zsrdg2 - srdg2(ji,jj) ) * rhoic * r1_rdtice 1130 1131 !rdm_ice(ji,jj) = rdm_ice(ji,jj) + vsw(ji,jj) * rhoic ! gurvan: increase in ice volume du to seawater frozen in voids 1002 sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ji,jj) * rhoic * r1_rdtice 1003 wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ji,jj) * rhoic * r1_rdtice ! increase in ice volume du to seawater frozen in voids 1132 1004 1133 1005 !------------------------------------ … … 1155 1027 ! ij looping 1-icells 1156 1028 1157 msnow_mlt(ji,jj) = msnow_mlt(ji,jj) + rhosn*vsrdg(ji,jj)*(1.0-fsnowrdg) & ! rafting included 1158 & + rhosn*vsrft(ji,jj)*(1.0-fsnowrft) 1159 1160 esnow_mlt(ji,jj) = esnow_mlt(ji,jj) + esrdg(ji,jj)*(1.0-fsnowrdg) & !rafting included 1161 & + esrft(ji,jj)*(1.0-fsnowrft) 1029 msnow_mlt(ji,jj) = msnow_mlt(ji,jj) + rhosn*vsrdg(ji,jj)*(1.0-rn_fsnowrdg) & ! rafting included 1030 & + rhosn*vsrft(ji,jj)*(1.0-rn_fsnowrft) 1031 1032 ! in J/m2 (same as e_s) 1033 esnow_mlt(ji,jj) = esnow_mlt(ji,jj) - esrdg(ji,jj)*(1.0-rn_fsnowrdg) & !rafting included 1034 & - esrft(ji,jj)*(1.0-rn_fsnowrft) 1162 1035 1163 1036 !----------------------------------------------------------------- … … 1172 1045 dhr2(ji,jj) = hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1) 1173 1046 1174 END DO ! ij1047 END DO 1175 1048 1176 1049 !-------------------------------------------------------------------- … … 1179 1052 !-------------------------------------------------------------------- 1180 1053 DO jk = 1, nlay_i 1181 !CDIR NODEP1182 1054 DO ij = 1, icells 1183 1055 ji = indxi(ij) 1184 1056 jj = indxj(ij) 1185 1057 ! heat content of ridged ice 1186 erdg1(ji,jj,jk) = eicen_init(ji,jj,jk,jl1) * afrac(ji,jj) / ( 1._wp + ridge_por )1058 erdg1(ji,jj,jk) = eicen_init(ji,jj,jk,jl1) * afrac(ji,jj) 1187 1059 eirft(ji,jj,jk) = eicen_init(ji,jj,jk,jl1) * afrft(ji,jj) 1188 1060 e_i (ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - erdg1(ji,jj,jk) - eirft(ji,jj,jk) 1189 ! sea water heat content 1190 ztmelts = - tmut * sss_m(ji,jj) + rtt 1191 ! heat content per unit volume 1192 zdummy0 = - rcp * ( sst_m(ji,jj) + rt0 - rtt ) * vsw(ji,jj) 1193 1194 ! corrected sea water salinity 1195 zindb = MAX( 0._wp , SIGN( 1._wp , vsw(ji,jj) - epsi20 ) ) 1196 zdummy = zindb * ( srdg1(ji,jj) - srdg2(ji,jj) ) / MAX( ridge_por * vsw(ji,jj), epsi20 ) 1197 1198 ztmelts = - tmut * zdummy + rtt 1199 ersw(ji,jj,jk) = - rcp * ( ztmelts - rtt ) * vsw(ji,jj) 1200 1201 ! heat flux 1202 fheat_mec(ji,jj) = fheat_mec(ji,jj) + ( ersw(ji,jj,jk) - zdummy0 ) * r1_rdtice 1203 1204 ! Correct dimensions to avoid big values 1205 ersw(ji,jj,jk) = ersw(ji,jj,jk) * 1.e-09 1206 1207 ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 1208 ersw (ji,jj,jk) = ersw(ji,jj,jk) * area(ji,jj) * vsw(ji,jj) / REAL( nlay_i ) 1209 1061 1062 1063 ! enthalpy of the trapped seawater (J/m2, >0) 1064 ! clem: if sst>0, then ersw <0 (is that possible?) 1065 ersw(ji,jj,jk) = - rhoic * vsw(ji,jj) * rcp * sst_m(ji,jj) * r1_nlay_i 1066 1067 ! heat flux to the ocean 1068 hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ji,jj,jk) * r1_rdtice ! > 0 [W.m-2] ocean->ice flux 1069 1070 ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean 1210 1071 erdg2(ji,jj,jk) = erdg1(ji,jj,jk) + ersw(ji,jj,jk) 1211 END DO ! ij 1212 END DO !jk 1072 1073 END DO 1074 END DO 1213 1075 1214 1076 1215 1077 IF( con_i ) THEN 1216 1078 DO jk = 1, nlay_i 1217 !CDIR NODEP1218 1079 DO ij = 1, icells 1219 1080 ji = indxi(ij) 1220 1081 jj = indxj(ij) 1221 1082 eice_init(ji,jj) = eice_init(ji,jj) + erdg2(ji,jj,jk) - erdg1(ji,jj,jk) 1222 END DO ! ij 1223 END DO !jk 1224 ENDIF 1225 1226 IF( large_afrac ) THEN ! there is a bug 1227 !CDIR NODEP 1228 DO ij = 1, icells 1229 ji = indxi(ij) 1230 jj = indxj(ij) 1231 IF( afrac(ji,jj) > kamax + epsi10 ) THEN 1232 WRITE(numout,*) '' 1233 WRITE(numout,*) ' ardg > a_i' 1234 WRITE(numout,*) ' ardg, aicen_init : ', ardg1(ji,jj), aicen_init(ji,jj,jl1) 1235 ENDIF 1236 END DO 1237 ENDIF 1238 IF( large_afrft ) THEN ! there is a bug 1239 !CDIR NODEP 1240 DO ij = 1, icells 1241 ji = indxi(ij) 1242 jj = indxj(ij) 1243 IF( afrft(ji,jj) > kamax + epsi10 ) THEN 1244 WRITE(numout,*) '' 1245 WRITE(numout,*) ' arft > a_i' 1246 WRITE(numout,*) ' arft, aicen_init : ', arft1(ji,jj), aicen_init(ji,jj,jl1) 1247 ENDIF 1083 END DO 1248 1084 END DO 1249 1085 ENDIF … … 1253 1089 !------------------------------------------------------------------------------- 1254 1090 ! jl1 looping 1-jpl 1255 DO jl2 = ice_cat_bounds(1,1), ice_cat_bounds(1,2)1091 DO jl2 = 1, jpl 1256 1092 ! over categories to which ridged ice is transferred 1257 !CDIR NODEP1258 1093 DO ij = 1, icells 1259 1094 ji = indxi(ij) … … 1264 1099 ! Transfer area, volume, and energy accordingly. 1265 1100 1266 IF( hrmin(ji,jj,jl1) >= hi_max(jl2) .OR. & 1267 hrmax(ji,jj,jl1) <= hi_max(jl2-1) ) THEN 1101 IF( hrmin(ji,jj,jl1) >= hi_max(jl2) .OR. hrmax(ji,jj,jl1) <= hi_max(jl2-1) ) THEN 1268 1102 hL = 0._wp 1269 1103 hR = 0._wp … … 1279 1113 a_i (ji,jj ,jl2) = a_i (ji,jj ,jl2) + ardg2 (ji,jj) * farea 1280 1114 v_i (ji,jj ,jl2) = v_i (ji,jj ,jl2) + vrdg2 (ji,jj) * fvol(ji,jj) 1281 v_s (ji,jj ,jl2) = v_s (ji,jj ,jl2) + vsrdg (ji,jj) * fvol(ji,jj) * fsnowrdg1282 e_s (ji,jj,1,jl2) = e_s (ji,jj,1,jl2) + esrdg (ji,jj) * fvol(ji,jj) * fsnowrdg1115 v_s (ji,jj ,jl2) = v_s (ji,jj ,jl2) + vsrdg (ji,jj) * fvol(ji,jj) * rn_fsnowrdg 1116 e_s (ji,jj,1,jl2) = e_s (ji,jj,1,jl2) + esrdg (ji,jj) * fvol(ji,jj) * rn_fsnowrdg 1283 1117 smv_i(ji,jj ,jl2) = smv_i(ji,jj ,jl2) + srdg2 (ji,jj) * fvol(ji,jj) 1284 1118 oa_i (ji,jj ,jl2) = oa_i (ji,jj ,jl2) + oirdg2(ji,jj) * farea 1285 1119 1286 END DO ! ij1120 END DO 1287 1121 1288 1122 ! Transfer ice energy to category jl2 by ridging 1289 1123 DO jk = 1, nlay_i 1290 !CDIR NODEP1291 1124 DO ij = 1, icells 1292 1125 ji = indxi(ij) 1293 1126 jj = indxj(ij) 1294 e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + fvol(ji,jj) *erdg2(ji,jj,jk)1127 e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + fvol(ji,jj) * erdg2(ji,jj,jk) 1295 1128 END DO 1296 1129 END DO … … 1298 1131 END DO ! jl2 (new ridges) 1299 1132 1300 DO jl2 = ice_cat_bounds(1,1), ice_cat_bounds(1,2) 1301 1302 !CDIR NODEP 1133 DO jl2 = 1, jpl 1134 1303 1135 DO ij = 1, icells 1304 1136 ji = indxi(ij) … … 1307 1139 ! thickness category jl2, transfer area, volume, and energy accordingly. 1308 1140 ! 1309 IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. & 1310 hraft(ji,jj,jl1) > hi_max(jl2-1) ) THEN 1141 IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) > hi_max(jl2-1) ) THEN 1311 1142 a_i (ji,jj ,jl2) = a_i (ji,jj ,jl2) + arft2 (ji,jj) 1312 1143 v_i (ji,jj ,jl2) = v_i (ji,jj ,jl2) + virft (ji,jj) 1313 v_s (ji,jj ,jl2) = v_s (ji,jj ,jl2) + vsrft (ji,jj) * fsnowrft1314 e_s (ji,jj,1,jl2) = e_s (ji,jj,1,jl2) + esrft (ji,jj) * fsnowrft1144 v_s (ji,jj ,jl2) = v_s (ji,jj ,jl2) + vsrft (ji,jj) * rn_fsnowrft 1145 e_s (ji,jj,1,jl2) = e_s (ji,jj,1,jl2) + esrft (ji,jj) * rn_fsnowrft 1315 1146 smv_i(ji,jj ,jl2) = smv_i(ji,jj ,jl2) + smrft (ji,jj) 1316 oa_i (ji,jj ,jl2) = oa_i (ji,jj ,jl2) + oirft2(ji,jj) 1317 ENDIF ! hraft1147 oa_i (ji,jj ,jl2) = oa_i (ji,jj ,jl2) + oirft2(ji,jj) 1148 ENDIF 1318 1149 ! 1319 END DO ! ij1150 END DO 1320 1151 1321 1152 ! Transfer rafted ice energy to category jl2 1322 1153 DO jk = 1, nlay_i 1323 !CDIR NODEP1324 1154 DO ij = 1, icells 1325 1155 ji = indxi(ij) 1326 1156 jj = indxj(ij) 1327 IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. & 1328 hraft(ji,jj,jl1) > hi_max(jl2-1) ) THEN 1157 IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) > hi_max(jl2-1) ) THEN 1329 1158 e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + eirft(ji,jj,jk) 1330 1159 ENDIF 1331 END DO ! ij1332 END DO !jk1333 1334 END DO ! jl21160 END DO 1161 END DO 1162 1163 END DO 1335 1164 1336 1165 END DO ! jl1 (deforming categories) … … 1346 1175 CALL lim_cons_check (eice_init, eice_final, 1.0e-2, fieldid) 1347 1176 1348 DO ji = mi0( jiindx), mi1(jiindx)1349 DO jj = mj0(j jindx), mj1(jjindx)1177 DO ji = mi0(iiceprt), mi1(iiceprt) 1178 DO jj = mj0(jiceprt), mj1(jiceprt) 1350 1179 WRITE(numout,*) ' vice_init : ', vice_init (ji,jj) 1351 1180 WRITE(numout,*) ' vice_final : ', vice_final(ji,jj) … … 1356 1185 ENDIF 1357 1186 ! 1358 CALL wrk_dealloc( (jpi+1)*(jpj+1), indxi, indxj )1359 CALL wrk_dealloc( jpi, jpj, vice_init, vice_final, eice_init, eice_final )1360 CALL wrk_dealloc( jpi, jpj, afrac, fvol , ardg1, ardg2, vsrdg, esrdg, oirdg1, oirdg2, dhr, dhr2 )1361 CALL wrk_dealloc( jpi, jpj, vrdg1, vrdg2, vsw , srdg1, srdg2, smsw)1362 CALL wrk_dealloc( jpi, jpj, afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )1363 CALL wrk_dealloc( jpi, jpj, jpl, aicen_init, vicen_init, vsnon_init, esnon_init, smv_i_init, oa_i_init )1364 CALL wrk_dealloc( jpi, jpj, jkmax,eirft, erdg1, erdg2, ersw )1365 CALL wrk_dealloc( jpi, jpj, jkmax, jpl,eicen_init )1187 CALL wrk_dealloc( (jpi+1)*(jpj+1), indxi, indxj ) 1188 CALL wrk_dealloc( jpi, jpj, vice_init, vice_final, eice_init, eice_final ) 1189 CALL wrk_dealloc( jpi, jpj, afrac, fvol , ardg1, ardg2, vsrdg, esrdg, dhr, dhr2 ) 1190 CALL wrk_dealloc( jpi, jpj, vrdg1, vrdg2, vsw , srdg1, srdg2, smsw, oirdg1, oirdg2 ) 1191 CALL wrk_dealloc( jpi, jpj, afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 ) 1192 CALL wrk_dealloc( jpi, jpj, jpl, aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init ) 1193 CALL wrk_dealloc( jpi, jpj, nlay_i, eirft, erdg1, erdg2, ersw ) 1194 CALL wrk_dealloc( jpi, jpj, nlay_i, jpl, eicen_init ) 1366 1195 ! 1367 1196 END SUBROUTINE lim_itd_me_ridgeshift 1368 1369 1370 SUBROUTINE lim_itd_me_asumr1371 !!-----------------------------------------------------------------------------1372 !! *** ROUTINE lim_itd_me_asumr ***1373 !!1374 !! ** Purpose : finds total fractional area1375 !!1376 !! ** Method : Find the total area of ice plus open water in each grid cell.1377 !! This is similar to the aggregate_area subroutine except that the1378 !! total area can be greater than 1, so the open water area is1379 !! included in the sum instead of being computed as a residual.1380 !!-----------------------------------------------------------------------------1381 INTEGER :: jl ! dummy loop index1382 !!-----------------------------------------------------------------------------1383 !1384 asum(:,:) = ato_i(:,:) ! open water1385 DO jl = 1, jpl ! ice categories1386 asum(:,:) = asum(:,:) + a_i(:,:,jl)1387 END DO1388 !1389 END SUBROUTINE lim_itd_me_asumr1390 1391 1197 1392 1198 SUBROUTINE lim_itd_me_init … … 1404 1210 !!------------------------------------------------------------------- 1405 1211 INTEGER :: ios ! Local integer output status for namelist read 1406 NAMELIST/namiceitdme/ ridge_scheme_swi, Cs, Cf, fsnowrdg, fsnowrft,& 1407 Gstar, astar, & 1408 Hstar, raftswi, hparmeter, Craft, ridge_por, & 1409 sal_max_ridge, partfun_swi, transfun_swi, & 1410 brinstren_swi 1212 NAMELIST/namiceitdme/ rn_cs, rn_fsnowrdg, rn_fsnowrft, & 1213 & rn_gstar, rn_astar, rn_hstar, ln_rafting, rn_hraft, rn_craft, rn_por_rdg, & 1214 & nn_partfun 1411 1215 !!------------------------------------------------------------------- 1412 1216 ! … … 1424 1228 WRITE(numout,*)' lim_itd_me_init : ice parameters for mechanical ice redistribution ' 1425 1229 WRITE(numout,*)' ~~~~~~~~~~~~~~~' 1426 WRITE(numout,*)' Switch choosing the ice redistribution scheme ridge_scheme_swi', ridge_scheme_swi 1427 WRITE(numout,*)' Fraction of shear energy contributing to ridging Cs ', Cs 1428 WRITE(numout,*)' Ratio of ridging work to PotEner change in ridging Cf ', Cf 1429 WRITE(numout,*)' Fraction of snow volume conserved during ridging fsnowrdg ', fsnowrdg 1430 WRITE(numout,*)' Fraction of snow volume conserved during ridging fsnowrft ', fsnowrft 1431 WRITE(numout,*)' Fraction of total ice coverage contributing to ridging Gstar ', Gstar 1432 WRITE(numout,*)' Equivalent to G* for an exponential part function astar ', astar 1433 WRITE(numout,*)' Quantity playing a role in max ridged ice thickness Hstar ', Hstar 1434 WRITE(numout,*)' Rafting of ice sheets or not raftswi ', raftswi 1435 WRITE(numout,*)' Parmeter thickness (threshold between ridge-raft) hparmeter ', hparmeter 1436 WRITE(numout,*)' Rafting hyperbolic tangent coefficient Craft ', Craft 1437 WRITE(numout,*)' Initial porosity of ridges ridge_por ', ridge_por 1438 WRITE(numout,*)' Maximum salinity of ridging ice sal_max_ridge ', sal_max_ridge 1439 WRITE(numout,*)' Switch for part. function (0) linear (1) exponential partfun_swi ', partfun_swi 1440 WRITE(numout,*)' Switch for tran. function (0) linear (1) exponential transfun_swi ', transfun_swi 1441 WRITE(numout,*)' Switch for including brine volume in ice strength comp. brinstren_swi ', brinstren_swi 1230 WRITE(numout,*)' Fraction of shear energy contributing to ridging rn_cs = ', rn_cs 1231 WRITE(numout,*)' Fraction of snow volume conserved during ridging rn_fsnowrdg = ', rn_fsnowrdg 1232 WRITE(numout,*)' Fraction of snow volume conserved during ridging rn_fsnowrft = ', rn_fsnowrft 1233 WRITE(numout,*)' Fraction of total ice coverage contributing to ridging rn_gstar = ', rn_gstar 1234 WRITE(numout,*)' Equivalent to G* for an exponential part function rn_astar = ', rn_astar 1235 WRITE(numout,*)' Quantity playing a role in max ridged ice thickness rn_hstar = ', rn_hstar 1236 WRITE(numout,*)' Rafting of ice sheets or not ln_rafting = ', ln_rafting 1237 WRITE(numout,*)' Parmeter thickness (threshold between ridge-raft) rn_hraft = ', rn_hraft 1238 WRITE(numout,*)' Rafting hyperbolic tangent coefficient rn_craft = ', rn_craft 1239 WRITE(numout,*)' Initial porosity of ridges rn_por_rdg = ', rn_por_rdg 1240 WRITE(numout,*)' Switch for part. function (0) linear (1) exponential nn_partfun = ', nn_partfun 1442 1241 ENDIF 1443 1242 ! 1444 1243 END SUBROUTINE lim_itd_me_init 1445 1446 1447 SUBROUTINE lim_itd_me_zapsmall1448 !!-------------------------------------------------------------------1449 !! *** ROUTINE lim_itd_me_zapsmall ***1450 !!1451 !! ** Purpose : Remove too small sea ice areas and correct salt fluxes1452 !!1453 !! history :1454 !! author: William H. Lipscomb, LANL1455 !! Nov 2003: Modified by Julie Schramm to conserve volume and energy1456 !! Sept 2004: Modified by William Lipscomb; replaced normalize_state with1457 !! additions to local freshwater, salt, and heat fluxes1458 !! 9.0, LIM3.0 - 02-2006 (M. Vancoppenolle) original code1459 !!-------------------------------------------------------------------1460 INTEGER :: ji, jj, jl, jk ! dummy loop indices1461 INTEGER :: icells ! number of cells with ice to zap1462 1463 REAL(wp), POINTER, DIMENSION(:,:) :: zmask ! 2D workspace1464 REAL(wp) :: zmask_glo1465 !!gm REAL(wp) :: xtmp ! temporary variable1466 !!-------------------------------------------------------------------1467 1468 CALL wrk_alloc( jpi, jpj, zmask )1469 1470 DO jl = 1, jpl1471 1472 !-----------------------------------------------------------------1473 ! Count categories to be zapped.1474 ! Abort model in case of negative area.1475 !-----------------------------------------------------------------1476 icells = 01477 zmask(:,:) = 0._wp1478 DO jj = 1, jpj1479 DO ji = 1, jpi1480 IF( ( a_i(ji,jj,jl) >= -epsi10 .AND. a_i(ji,jj,jl) < 0._wp ) .OR. &1481 & ( a_i(ji,jj,jl) > 0._wp .AND. a_i(ji,jj,jl) <= epsi10 ) .OR. &1482 & ( v_i(ji,jj,jl) == 0._wp .AND. a_i(ji,jj,jl) > 0._wp ) .OR. &1483 & ( v_i(ji,jj,jl) > 0._wp .AND. v_i(ji,jj,jl) <= epsi10 ) ) zmask(ji,jj) = 1._wp1484 END DO1485 END DO1486 !zmask_glo = glob_sum(zmask)1487 !IF( ln_nicep .AND. lwp ) WRITE(numout,*) zmask_glo, ' cells of ice zapped in the ocean '1488 1489 !-----------------------------------------------------------------1490 ! Zap ice energy and use ocean heat to melt ice1491 !-----------------------------------------------------------------1492 1493 DO jk = 1, nlay_i1494 DO jj = 1 , jpj1495 DO ji = 1 , jpi1496 !!gm xtmp = e_i(ji,jj,jk,jl) / area(ji,jj) * r1_rdtice1497 !!gm xtmp = xtmp * unit_fac1498 ! fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp1499 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * ( 1._wp - zmask(ji,jj) )1500 END DO1501 END DO1502 END DO1503 1504 DO jj = 1 , jpj1505 DO ji = 1 , jpi1506 1507 !-----------------------------------------------------------------1508 ! Zap snow energy and use ocean heat to melt snow1509 !-----------------------------------------------------------------1510 ! xtmp = esnon(i,j,n) / dt ! < 01511 ! fhnet(i,j) = fhnet(i,j) + xtmp1512 ! fhnet_hist(i,j) = fhnet_hist(i,j) + xtmp1513 ! xtmp is greater than 01514 ! fluxes are positive to the ocean1515 ! here the flux has to be negative for the ocean1516 !!gm xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) * r1_rdtice1517 ! fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp1518 1519 !!gm xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) * r1_rdtice !RB ???????1520 1521 t_s(ji,jj,1,jl) = rtt * zmask(ji,jj) + t_s(ji,jj,1,jl) * ( 1._wp - zmask(ji,jj) )1522 1523 !-----------------------------------------------------------------1524 ! zap ice and snow volume, add water and salt to ocean1525 !-----------------------------------------------------------------1526 1527 ! xtmp = (rhoi*vicen(i,j,n) + rhos*vsnon(i,j,n)) / dt1528 ! sfx_res(ji,jj) = sfx_res(ji,jj) + ( sss_m(ji,jj) ) &1529 ! * rhosn * v_s(ji,jj,jl) * r1_rdtice1530 ! sfx_res(ji,jj) = sfx_res(ji,jj) + ( sss_m(ji,jj) - sm_i(ji,jj,jl) ) &1531 ! * rhoic * v_i(ji,jj,jl) * r1_rdtice1532 ! sfx (i,j) = sfx (i,j) + xtmp1533 1534 ato_i(ji,jj) = a_i (ji,jj,jl) * zmask(ji,jj) + ato_i(ji,jj)1535 a_i (ji,jj,jl) = a_i (ji,jj,jl) * ( 1._wp - zmask(ji,jj) )1536 v_i (ji,jj,jl) = v_i (ji,jj,jl) * ( 1._wp - zmask(ji,jj) )1537 v_s (ji,jj,jl) = v_s (ji,jj,jl) * ( 1._wp - zmask(ji,jj) )1538 t_su (ji,jj,jl) = t_su (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) + t_bo(ji,jj) * zmask(ji,jj)1539 oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * ( 1._wp - zmask(ji,jj) )1540 smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * ( 1._wp - zmask(ji,jj) )1541 !1542 END DO1543 END DO1544 !1545 END DO ! jl1546 !1547 CALL wrk_dealloc( jpi, jpj, zmask )1548 !1549 END SUBROUTINE lim_itd_me_zapsmall1550 1244 1551 1245 #else … … 1558 1252 SUBROUTINE lim_itd_me_icestrength 1559 1253 END SUBROUTINE lim_itd_me_icestrength 1560 SUBROUTINE lim_itd_me_sort1561 END SUBROUTINE lim_itd_me_sort1562 1254 SUBROUTINE lim_itd_me_init 1563 1255 END SUBROUTINE lim_itd_me_init 1564 SUBROUTINE lim_itd_me_zapsmall1565 END SUBROUTINE lim_itd_me_zapsmall1566 1256 #endif 1567 1257 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
