Changeset 990 for branches/dev_003_CPL/NEMO/LIM_SRC_3/limthd_dh.F90
- Timestamp:
- 2008-05-23T16:38:21+02:00 (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/dev_003_CPL/NEMO/LIM_SRC_3/limthd_dh.F90
r888 r990 24 24 USE par_ice 25 25 USE lib_mpp 26 26 27 27 IMPLICIT NONE 28 28 PRIVATE … … 46 46 47 47 SUBROUTINE lim_thd_dh(kideb,kiut,jl) 48 !!------------------------------------------------------------------ 49 !! *** ROUTINE lim_thd_dh *** 50 !!------------------------------------------------------------------ 51 !! ** Purpose : 52 !! This routine determines variations of ice and snow thicknesses. 53 !! ** Method : 54 !! Ice/Snow surface melting arises from imbalance in surface fluxes 55 !! Bottom accretion/ablation arises from flux budget 56 !! Snow thickness can increase by precipitation and decrease by 57 !! sublimation 58 !! If snow load excesses Archmiede limit, snow-ice is formed by 59 !! the flooding of sea-water in the snow 60 !! ** Steps 61 !! 1) Compute available flux of heat for surface ablation 62 !! 2) Compute snow and sea ice enthalpies 63 !! 3) Surface ablation and sublimation 64 !! 4) Bottom accretion/ablation 65 !! 5) Case of Total ablation 66 !! 6) Snow ice formation 67 !! 68 !! ** Arguments 69 !! 70 !! ** Inputs / Outputs 71 !! 72 !! ** External 73 !! 74 !! ** References : Bitz and Lipscomb, JGR 99 75 !! Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646 76 !! Vancoppenolle, Fichefet and Bitz, GRL 2005 77 !! Vancoppenolle et al., OM08 78 !! 79 !! ** History : 80 !! original code 01-04 (LIM) 81 !! original routine 82 !! (05-2003) M. Vancoppenolle, Louvain-La-Neuve, Belgium 83 !! (06/07-2005) 3D version 84 !! (03-2008) Clean code 85 !! 86 !!------------------------------------------------------------------ 87 !! * Arguments 88 INTEGER , INTENT (IN) :: & 89 kideb , & !: Start point on which the the computation is applied 90 kiut , & !: End point on which the the computation is applied 91 jl !: Thickness cateogry number 92 93 !! * Local variables 94 INTEGER :: & 95 ji , & !: space index 96 jk , & !: ice layer index 97 isnow , & !: switch for presence (1) or absence (0) of snow 98 zji , & !: 2D corresponding indices to ji 99 zjj , & 100 isnowic , & !: snow ice formation not 101 i_ice_switch , & !: ice thickness above a certain treshold or not 102 iter 103 104 REAL(wp) :: & 105 zhsnew , & !: new snow thickness 106 zihgnew , & !: switch for total ablation 107 ztmelts , & !: melting point 108 zhn , & 109 zdhcf , & 110 zdhbf , & 111 zhni , & 112 zhnfi , & 113 zihg , & 114 zdhnm , & 115 zhnnew , & 116 zeps = 1.0e-13, & 117 zhisn , & 118 zfracs , & !: fractionation coefficient for bottom salt 119 !: entrapment 120 zds , & !: increment of bottom ice salinity 121 zcoeff , & !: dummy argument for snowfall partitioning 122 !: over ice and leads 123 zsm_snowice, & !: snow-ice salinity 124 zswi1 , & !: switch for computation of bottom salinity 125 zswi12 , & !: switch for computation of bottom salinity 126 zswi2 , & !: switch for computation of bottom salinity 127 zgrr , & !: bottom growth rate 128 zihic , & !: 129 ztform !: bottom formation temperature 130 131 REAL(wp) , DIMENSION(jpij) :: & 132 zh_i , & ! ice layer thickness 133 zh_s , & ! snow layer thickness 134 ztfs , & ! melting point 135 zhsold , & ! old snow thickness 136 zqprec , & !: energy of fallen snow 137 zqfont_su , & ! incoming, remaining surface energy 138 zqfont_bo ! incoming, bottom energy 139 140 REAL(wp) , DIMENSION(jpij) :: & 141 z_f_surf, & ! surface heat for ablation 142 zhgnew ! new ice thickness 143 144 REAL(wp), DIMENSION(jpij) :: & 145 zdh_s_mel , & ! snow melt 146 zdh_s_pre , & ! snow precipitation 147 zdh_s_sub , & ! snow sublimation 148 zfsalt_melt ! salt flux due to ice melt 149 150 REAL(wp) , DIMENSION(jpij,jkmax) :: & 151 zdeltah 152 153 ! Pathological cases 154 REAL(wp), DIMENSION(jpij) :: & 155 zfdt_init , & !: total incoming heat for ice melt 156 zfdt_final , & !: total remaing heat for ice melt 157 zqt_i , & !: total ice heat content 158 zqt_s , & !: total snow heat content 159 zqt_dummy !: dummy heat content 160 161 REAL(wp), DIMENSION(jpij,jkmax) :: & 162 zqt_i_lay !: total ice heat content 163 164 ! Heat conservation 165 REAL(wp), DIMENSION(jpij) :: & 166 zfbase, & 167 zdq_i 168 169 INTEGER, DIMENSION(jpij) :: & 170 innermelt 171 172 REAL(wp) :: & 173 meance_dh 174 175 INTEGER :: & 176 num_iter_max, & 177 numce_dh 178 179 !!----------------------------------------------------------------------------- 180 181 WRITE(numout,*) 'lim_thd_dh : computation of vertical snow/ice accretion/ablation' 182 WRITE(numout,*) '~~~~~~~~~' 48 !!------------------------------------------------------------------ 49 !! *** ROUTINE lim_thd_dh *** 50 !!------------------------------------------------------------------ 51 !! ** Purpose : 52 !! This routine determines variations of ice and snow thicknesses. 53 !! ** Method : 54 !! Ice/Snow surface melting arises from imbalance in surface fluxes 55 !! Bottom accretion/ablation arises from flux budget 56 !! Snow thickness can increase by precipitation and decrease by 57 !! sublimation 58 !! If snow load excesses Archmiede limit, snow-ice is formed by 59 !! the flooding of sea-water in the snow 60 !! ** Steps 61 !! 1) Compute available flux of heat for surface ablation 62 !! 2) Compute snow and sea ice enthalpies 63 !! 3) Surface ablation and sublimation 64 !! 4) Bottom accretion/ablation 65 !! 5) Case of Total ablation 66 !! 6) Snow ice formation 67 !! 68 !! ** Arguments 69 !! 70 !! ** Inputs / Outputs 71 !! 72 !! ** External 73 !! 74 !! ** References : Bitz and Lipscomb, JGR 99 75 !! Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646 76 !! Vancoppenolle, Fichefet and Bitz, GRL 2005 77 !! Vancoppenolle et al., OM08 78 !! 79 !! ** History : 80 !! original code 01-04 (LIM) 81 !! original routine 82 !! (05-2003) M. Vancoppenolle, Louvain-La-Neuve, Belgium 83 !! (06/07-2005) 3D version 84 !! (03-2008) Clean code 85 !! 86 !!------------------------------------------------------------------ 87 !! * Arguments 88 INTEGER , INTENT (IN) :: & 89 kideb , & !: Start point on which the the computation is applied 90 kiut , & !: End point on which the the computation is applied 91 jl !: Thickness cateogry number 92 93 !! * Local variables 94 INTEGER :: & 95 ji , & !: space index 96 jk , & !: ice layer index 97 isnow , & !: switch for presence (1) or absence (0) of snow 98 zji , & !: 2D corresponding indices to ji 99 zjj , & 100 isnowic , & !: snow ice formation not 101 i_ice_switch , & !: ice thickness above a certain treshold or not 102 iter 103 104 REAL(wp) :: & 105 zhsnew , & !: new snow thickness 106 zihgnew , & !: switch for total ablation 107 ztmelts , & !: melting point 108 zhn , & 109 zdhcf , & 110 zdhbf , & 111 zhni , & 112 zhnfi , & 113 zihg , & 114 zdhnm , & 115 zhnnew , & 116 zeps = 1.0e-13, & 117 zhisn , & 118 zfracs , & !: fractionation coefficient for bottom salt 119 !: entrapment 120 zds , & !: increment of bottom ice salinity 121 zcoeff , & !: dummy argument for snowfall partitioning 122 !: over ice and leads 123 zsm_snowice, & !: snow-ice salinity 124 zswi1 , & !: switch for computation of bottom salinity 125 zswi12 , & !: switch for computation of bottom salinity 126 zswi2 , & !: switch for computation of bottom salinity 127 zgrr , & !: bottom growth rate 128 zihic , & !: 129 ztform !: bottom formation temperature 130 131 REAL(wp) , DIMENSION(jpij) :: & 132 zh_i , & ! ice layer thickness 133 zh_s , & ! snow layer thickness 134 ztfs , & ! melting point 135 zhsold , & ! old snow thickness 136 zqprec , & !: energy of fallen snow 137 zqfont_su , & ! incoming, remaining surface energy 138 zqfont_bo ! incoming, bottom energy 139 140 REAL(wp) , DIMENSION(jpij) :: & 141 z_f_surf, & ! surface heat for ablation 142 zhgnew ! new ice thickness 143 144 REAL(wp), DIMENSION(jpij) :: & 145 zdh_s_mel , & ! snow melt 146 zdh_s_pre , & ! snow precipitation 147 zdh_s_sub , & ! snow sublimation 148 zfsalt_melt ! salt flux due to ice melt 149 150 REAL(wp) , DIMENSION(jpij,jkmax) :: & 151 zdeltah 152 153 ! Pathological cases 154 REAL(wp), DIMENSION(jpij) :: & 155 zfdt_init , & !: total incoming heat for ice melt 156 zfdt_final , & !: total remaing heat for ice melt 157 zqt_i , & !: total ice heat content 158 zqt_s , & !: total snow heat content 159 zqt_dummy !: dummy heat content 160 161 REAL(wp), DIMENSION(jpij,jkmax) :: & 162 zqt_i_lay !: total ice heat content 163 164 ! Heat conservation 165 REAL(wp), DIMENSION(jpij) :: & 166 zfbase, & 167 zdq_i 168 169 INTEGER, DIMENSION(jpij) :: & 170 innermelt 171 172 REAL(wp) :: & 173 meance_dh 174 175 INTEGER :: & 176 num_iter_max, & 177 numce_dh 183 178 184 179 zfsalt_melt(:) = 0.0 … … 191 186 old_ht_s_b(ji) = ht_s_b(ji) 192 187 END DO 193 !194 !------------------------------------------------------------------------------!195 ! 1) Calculate available heat for surface ablation !196 !------------------------------------------------------------------------------!197 !188 ! 189 !------------------------------------------------------------------------------! 190 ! 1) Calculate available heat for surface ablation ! 191 !------------------------------------------------------------------------------! 192 ! 198 193 DO ji = kideb,kiut 199 194 isnow = INT( 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_s_b(ji) ) ) ) 200 195 ztfs(ji) = isnow * rtt + ( 1.0 - isnow ) * rtt 201 196 z_f_surf(ji) = qnsr_ice_1d(ji) + ( 1.0 - i0(ji) ) * & 202 197 qsr_ice_1d(ji) - fc_su(ji) 203 198 z_f_surf(ji) = MAX( zzero , z_f_surf(ji) ) * & 204 199 MAX( zzero , SIGN( zone , t_su_b(ji) - ztfs(ji) ) ) 205 200 zfdt_init(ji) = ( z_f_surf(ji) + & 206 207 201 MAX( fbif_1d(ji) + qlbbq_1d(ji) + fc_bo_i(ji),0.0 ) ) & 202 * rdt_ice 208 203 END DO ! ji 209 204 … … 212 207 dsm_i_se_1d(:) = 0.0 213 208 dsm_i_si_1d(:) = 0.0 214 !215 !------------------------------------------------------------------------------!216 ! 2) Computing layer thicknesses and snow and sea-ice enthalpies. !217 !------------------------------------------------------------------------------!218 !209 ! 210 !------------------------------------------------------------------------------! 211 ! 2) Computing layer thicknesses and snow and sea-ice enthalpies. ! 212 !------------------------------------------------------------------------------! 213 ! 219 214 ! Layer thickness 220 215 DO ji = kideb,kiut … … 239 234 END DO 240 235 END DO 241 !242 !------------------------------------------------------------------------------|243 ! 3) Surface ablation and sublimation |244 !------------------------------------------------------------------------------|245 !236 ! 237 !------------------------------------------------------------------------------| 238 ! 3) Surface ablation and sublimation | 239 !------------------------------------------------------------------------------| 240 ! 246 241 !------------------------- 247 242 ! 3.1 Snow precips / melt … … 272 267 zdeltah(ji,1) = MIN( 0.0 , - zqfont_su(ji) / MAX( zqprec(ji) , epsi13 ) ) 273 268 zqfont_su(ji) = MAX( 0.0 , - zdh_s_pre(ji) - zdeltah(ji,1) ) * & 274 269 zqprec(ji) 275 270 zdeltah(ji,1) = MAX( - zdh_s_pre(ji) , zdeltah(ji,1) ) 276 271 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,1) … … 289 284 zdeltah(ji,jk) = - zqfont_su(ji) / q_s_b(ji,jk) 290 285 zqfont_su(ji) = MAX( 0.0 , - zh_s(ji) - zdeltah(ji,jk) ) * & 291 286 q_s_b(ji,jk) 292 287 zdeltah(ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) ) 293 288 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) !resulting melt of snow … … 306 301 ! Volume and mass variations of snow 307 302 dvsbq_1d(ji) = a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji) & 308 303 - zdh_s_mel(ji) ) 309 304 dvsbq_1d(ji) = MIN( zzero, dvsbq_1d(ji) ) 310 305 rdmsnif_1d(ji) = rhosn*dvsbq_1d(ji) … … 327 322 ! recompute heat available 328 323 zqfont_su(ji) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * & 329 324 q_i_b(ji,jk) 330 325 ! melt of layer jk cannot be higher than its thickness 331 326 zdeltah(ji,jk) = MAX( zdeltah(ji,jk) , - zh_i(ji) ) … … 334 329 ! for energy conservation 335 330 zdq_i(ji) = zdq_i(ji) + zdeltah(ji,jk) * & 336 331 q_i_b(ji,jk) / rdt_ice 337 332 ! contribution to ice-ocean salt flux 338 333 zji = MOD( npb(ji) - 1, jpi ) + 1 339 334 zjj = ( npb(ji) - 1 ) / jpi + 1 340 335 zfsalt_melt(ji) = zfsalt_melt(ji) + & 341 342 343 336 ( sss_m(zji,zjj) - sm_i_b(ji) ) * & 337 a_i_b(ji) * & 338 MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice 344 339 END DO ! ji 345 340 END DO ! jk … … 349 344 !------------------- 350 345 IF ( con_i ) THEN 351 numce_dh = 0352 meance_dh = 0.0353 DO ji = kideb, kiut354 355 IF ( ( z_f_surf(ji) + zdq_i(ji) ) .GE. 1.0e-3 ) THEN356 numce_dh = numce_dh + 1357 meance_dh = meance_dh + z_f_surf(ji) + zdq_i(ji)358 ENDIF359 360 IF ( z_f_surf(ji) + zdq_i(ji) .GE. 1.0e-3 ) THEN!361 WRITE(numout,*) ' ALERTE heat loss for surface melt '362 WRITE(numout,*) ' zji, zjj, jl :', zji, zjj, jl363 WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji)364 WRITE(numout,*) ' z_f_surf : ', z_f_surf(ji)365 WRITE(numout,*) ' zdq_i : ', zdq_i(ji)366 WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji)367 WRITE(numout,*) ' fc_bo_i : ', fc_bo_i(ji)368 WRITE(numout,*) ' fbif_1d : ', fbif_1d(ji)369 WRITE(numout,*) ' qlbbq_1d: ', qlbbq_1d(ji)370 WRITE(numout,*) ' s_i_new : ', s_i_new(ji)371 WRITE(numout,*) ' sss_m : ', sss_m(zji,zjj)372 ENDIF373 374 END DO ! ji375 376 IF ( numce_dh .GT. 0 ) meance_dh = meance_dh / numce_dh377 WRITE(numout,*) ' Error report - Category : ', jl378 WRITE(numout,*) ' ~~~~~~~~~~~~ '379 WRITE(numout,*) ' Number of points where there is sur. me. error : ', numce_dh380 WRITE(numout,*) ' Mean basal growth error on error points : ', meance_dh346 numce_dh = 0 347 meance_dh = 0.0 348 DO ji = kideb, kiut 349 350 IF ( ( z_f_surf(ji) + zdq_i(ji) ) .GE. 1.0e-3 ) THEN 351 numce_dh = numce_dh + 1 352 meance_dh = meance_dh + z_f_surf(ji) + zdq_i(ji) 353 ENDIF 354 355 IF ( z_f_surf(ji) + zdq_i(ji) .GE. 1.0e-3 ) THEN! 356 WRITE(numout,*) ' ALERTE heat loss for surface melt ' 357 WRITE(numout,*) ' zji, zjj, jl :', zji, zjj, jl 358 WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji) 359 WRITE(numout,*) ' z_f_surf : ', z_f_surf(ji) 360 WRITE(numout,*) ' zdq_i : ', zdq_i(ji) 361 WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji) 362 WRITE(numout,*) ' fc_bo_i : ', fc_bo_i(ji) 363 WRITE(numout,*) ' fbif_1d : ', fbif_1d(ji) 364 WRITE(numout,*) ' qlbbq_1d: ', qlbbq_1d(ji) 365 WRITE(numout,*) ' s_i_new : ', s_i_new(ji) 366 WRITE(numout,*) ' sss_m : ', sss_m(zji,zjj) 367 ENDIF 368 369 END DO ! ji 370 371 IF ( numce_dh .GT. 0 ) meance_dh = meance_dh / numce_dh 372 WRITE(numout,*) ' Error report - Category : ', jl 373 WRITE(numout,*) ' ~~~~~~~~~~~~ ' 374 WRITE(numout,*) ' Number of points where there is sur. me. error : ', numce_dh 375 WRITE(numout,*) ' Mean basal growth error on error points : ', meance_dh 381 376 382 377 ENDIF ! con_i … … 409 404 DO jk = 1, nlay_s !n 410 405 DO ji = kideb, kiut !n 411 ! In case of disparition of the snow, we have to update the snow412 ! temperatures406 ! In case of disparition of the snow, we have to update the snow 407 ! temperatures 413 408 zhisn = MAX( zzero , SIGN( zone, - ht_s_b(ji) ) ) 414 409 t_s_b(ji,jk) = ( 1.0 - zhisn ) * t_s_b(ji,jk) + zhisn * rtt 415 410 q_s_b(ji,jk) = ( 1.0 - zhisn ) * q_s_b(ji,jk) 416 411 END DO 417 END DO 418 419 !420 !------------------------------------------------------------------------------!421 ! 4) Basal growth / melt !422 !------------------------------------------------------------------------------!423 !412 END DO 413 414 ! 415 !------------------------------------------------------------------------------! 416 ! 4) Basal growth / melt ! 417 !------------------------------------------------------------------------------! 418 ! 424 419 ! Ice basal growth / melt is given by the ratio of heat budget over basal 425 420 ! ice heat content. Basal heat budget is given by the difference between … … 439 434 ! New ice heat content (Bitz and Lipscomb, 1999) 440 435 ztform = t_i_b(ji,nlay_i) ! t_bo_b crashes in the 441 436 ! Baltic 442 437 q_i_b(ji,nlay_i+1) = rhoic * & 443 444 445 446 438 ( cpic * ( ztmelts - ztform ) & 439 + lfus * ( 1.0 - ( ztmelts - rtt ) / & 440 ( ztform - rtt ) ) & 441 - rcp * ( ztmelts-rtt ) ) 447 442 ! Basal growth rate = - F*dt / q 448 443 dh_i_bott(ji) = - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + & 449 444 qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 450 445 ENDIF ! heat budget 451 446 END DO ! ji … … 476 471 ! New ice heat content (Bitz and Lipscomb, 1999) 477 472 q_i_b(ji,nlay_i+1) = rhoic * & 478 479 480 481 473 ( cpic * ( ztmelts - t_bo_b(ji) ) & 474 + lfus * ( 1.0 - ( ztmelts - rtt ) / & 475 ( t_bo_b(ji) - rtt ) ) & 476 - rcp * ( ztmelts-rtt ) ) 482 477 ! Bottom growth rate = - F*dt / q 483 478 dh_i_bott(ji) = - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) & 484 479 + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 485 480 ! New ice salinity ( Cox and Weeks, JGR, 1988 ) 486 481 ! zswi2 (1) if dh_i_bott/rdt .GT. 3.6e-7 … … 492 487 zswi1 = 1. - zswi2 * zswi12 493 488 zfracs = zswi1 * 0.12 + & 494 495 496 489 zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) ) + & 490 zswi2 * 0.26 / & 491 ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) 497 492 zds = zfracs*sss_m(zji,zjj) - s_i_new(ji) 498 493 s_i_new(ji) = zfracs * sss_m(zji,zjj) … … 510 505 ! New ice heat content (Bitz and Lipscomb, 1999) 511 506 q_i_b(ji,nlay_i+1) = rhoic * & 512 513 514 515 507 ( cpic * ( ztmelts - t_bo_b(ji) ) & 508 + lfus * ( 1.0 - ( ztmelts - rtt ) / & 509 ( t_bo_b(ji) - rtt ) ) & 510 - rcp * ( ztmelts-rtt ) ) 516 511 ! Basal growth rate = - F*dt / q 517 512 dh_i_bott(ji) = - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + & 518 513 qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 519 514 ! Salinity update 520 515 ! entrapment during bottom growth 521 516 dsm_i_se_1d(ji) = ( s_i_new(ji)*dh_i_bott(ji) + & 522 523 524 517 sm_i_b(ji)*ht_i_b(ji) ) / & 518 MAX( ht_i_b(ji) + dh_i_bott(ji) ,zeps ) & 519 - sm_i_b(ji) 525 520 ENDIF ! heat budget 526 521 END DO ! ji … … 537 532 ! heat convergence at the surface > 0 538 533 IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0 ) THEN 539 534 540 535 s_i_new(ji) = s_i_b(ji,nlay_i) 541 536 zqfont_bo(ji) = rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) … … 559 554 zdeltah(ji,jk) = - zqfont_bo(ji) / q_i_b(ji,jk) 560 555 zqfont_bo(ji) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * & 561 556 q_i_b(ji,jk) 562 557 zdeltah(ji,jk) = MAX(zdeltah(ji,jk), - zh_i(ji) ) 563 558 dh_i_bott(ji) = dh_i_bott(ji) + zdeltah(ji,jk) 564 559 zdq_i(ji) = zdq_i(ji) + zdeltah(ji,jk) * & 565 566 ! contribution to salt flux560 q_i_b(ji,jk) / rdt_ice 561 ! contribution to salt flux 567 562 zji = MOD( npb(ji) - 1, jpi ) + 1 568 563 zjj = ( npb(ji) - 1 ) / jpi + 1 569 564 zfsalt_melt(ji) = zfsalt_melt(ji) + & 570 571 572 565 ( sss_m(zji,zjj) - sm_i_b(ji) ) * & 566 a_i_b(ji) * & 567 MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice 573 568 ENDIF 574 569 ENDIF … … 580 575 !------------------- 581 576 IF ( con_i ) THEN 582 DO ji = kideb, kiut583 IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0 ) THEN584 IF ( ( zfbase(ji) + zdq_i(ji) ) .GE. 1.0e-3 ) THEN585 numce_dh = numce_dh + 1586 meance_dh = meance_dh + zfbase(ji) + zdq_i(ji)587 ENDIF588 IF ( zfbase(ji) + zdq_i(ji) .GE. 1.0e-3 ) THEN589 WRITE(numout,*) ' ALERTE heat loss for basal melt '590 WRITE(numout,*) ' zji, zjj, jl :', zji, zjj, jl591 WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji)592 WRITE(numout,*) ' zfbase : ', zfbase(ji)593 WRITE(numout,*) ' zdq_i : ', zdq_i(ji)594 WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji)595 WRITE(numout,*) ' fc_bo_i : ', fc_bo_i(ji)596 WRITE(numout,*) ' fbif_1d : ', fbif_1d(ji)597 WRITE(numout,*) ' qlbbq_1d: ', qlbbq_1d(ji)598 WRITE(numout,*) ' s_i_new : ', s_i_new(ji)599 WRITE(numout,*) ' sss_m : ', sss_m(zji,zjj)600 WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)601 WRITE(numout,*) ' innermelt : ', innermelt(ji)602 ENDIF603 ENDIF ! heat convergence at the surface604 END DO ! ji605 606 IF ( numce_dh .GT. 0 ) meance_dh = meance_dh / numce_dh607 WRITE(numout,*) ' Number of points where there is bas. me. error : ', numce_dh608 WRITE(numout,*) ' Mean basal melt error on error points : ', meance_dh609 WRITE(numout,*) ' Remaining bottom heat : ', zqfont_bo(jiindex_1d)577 DO ji = kideb, kiut 578 IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0 ) THEN 579 IF ( ( zfbase(ji) + zdq_i(ji) ) .GE. 1.0e-3 ) THEN 580 numce_dh = numce_dh + 1 581 meance_dh = meance_dh + zfbase(ji) + zdq_i(ji) 582 ENDIF 583 IF ( zfbase(ji) + zdq_i(ji) .GE. 1.0e-3 ) THEN 584 WRITE(numout,*) ' ALERTE heat loss for basal melt ' 585 WRITE(numout,*) ' zji, zjj, jl :', zji, zjj, jl 586 WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji) 587 WRITE(numout,*) ' zfbase : ', zfbase(ji) 588 WRITE(numout,*) ' zdq_i : ', zdq_i(ji) 589 WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji) 590 WRITE(numout,*) ' fc_bo_i : ', fc_bo_i(ji) 591 WRITE(numout,*) ' fbif_1d : ', fbif_1d(ji) 592 WRITE(numout,*) ' qlbbq_1d: ', qlbbq_1d(ji) 593 WRITE(numout,*) ' s_i_new : ', s_i_new(ji) 594 WRITE(numout,*) ' sss_m : ', sss_m(zji,zjj) 595 WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji) 596 WRITE(numout,*) ' innermelt : ', innermelt(ji) 597 ENDIF 598 ENDIF ! heat convergence at the surface 599 END DO ! ji 600 601 IF ( numce_dh .GT. 0 ) meance_dh = meance_dh / numce_dh 602 WRITE(numout,*) ' Number of points where there is bas. me. error : ', numce_dh 603 WRITE(numout,*) ' Mean basal melt error on error points : ', meance_dh 604 WRITE(numout,*) ' Remaining bottom heat : ', zqfont_bo(jiindex_1d) 610 605 611 606 ENDIF ! con_i 612 607 613 !614 !------------------------------------------------------------------------------!615 ! 5) Pathological cases !616 !------------------------------------------------------------------------------!617 !608 ! 609 !------------------------------------------------------------------------------! 610 ! 5) Pathological cases ! 611 !------------------------------------------------------------------------------! 612 ! 618 613 !---------------------------------------------- 619 614 ! 5.1 Excessive ablation in a 1-category model … … 626 621 ! excessive energy is sent to lateral ablation 627 622 fsup(ji) = rhoic*lfus * at_i_b(ji) / MAX( ( 1.0 - at_i_b(ji) ),epsi13) & 628 623 * ( zdhbf - dh_i_bott(ji) ) / rdt_ice 629 624 630 625 dh_i_bott(ji) = zdhbf … … 638 633 zjj = ( npb(ji) - 1 ) / jpi + 1 639 634 diag_bot_gr(zji,zjj) = diag_bot_gr(zji,zjj) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) & 640 635 / rdt_ice 641 636 diag_sur_me(zji,zjj) = diag_sur_me(zji,zjj) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) & 642 637 / rdt_ice 643 638 diag_bot_me(zji,zjj) = diag_bot_me(zji,zjj) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) & 644 639 / rdt_ice 645 640 END DO 646 641 … … 667 662 zqt_s(ji) = ( 1. - zihg) * zqt_s(ji) / MAX( zhni, zeps ) 668 663 zdhnm = - ( 1. - zihg ) * ( 1. - zihgnew ) * ( zfdt_final(ji) / & 669 664 MAX( zqt_s(ji) , zeps ) ) 670 665 zhnfi = zhni + zdhnm 671 666 zfdt_final(ji) = MAX ( zfdt_final(ji) + zqt_s(ji) * zdhnm , 0.0 ) … … 676 671 !--------------------------------- 677 672 rdmicif_1d(ji) = rdmicif_1d(ji) + a_i_b(ji) * & 678 673 (zhgnew(ji)-ht_i_b(ji))*rhoic ! good 679 674 680 675 rdmsnif_1d(ji) = rdmsnif_1d(ji) + a_i_b(ji) * & 681 676 (ht_s_b(ji)-zhni)*rhosn ! good too 682 677 683 678 ! Remaining heat to the ocean … … 700 695 zjj = ( npb(ji) - 1 ) / jpi + 1 701 696 IF ( num_sal .NE. 4 ) & 702 fseqv_1d(ji) = fseqv_1d(ji) + zihgnew * zfsalt_melt(ji) + &703 704 697 fseqv_1d(ji) = fseqv_1d(ji) + zihgnew * zfsalt_melt(ji) + & 698 (1.0 - zihgnew) * rdmicif_1d(ji) * & 699 ( sss_m(zji,zjj) - sm_i_b(ji) ) / rdt_ice 705 700 ! new lines 706 701 IF ( num_sal .EQ. 4 ) & 707 fseqv_1d(ji) = fseqv_1d(ji) + zihgnew * zfsalt_melt(ji) + &708 709 702 fseqv_1d(ji) = fseqv_1d(ji) + zihgnew * zfsalt_melt(ji) + & 703 (1.0 - zihgnew) * rdmicif_1d(ji) * & 704 ( sss_m(zji,zjj) - bulk_sal ) / rdt_ice 710 705 ! Heat flux 711 706 ! excessive bottom ablation energy (fsup) - 0 except if jpl = 1 712 707 ! excessive total ablation energy (focea) sent to the ocean 713 708 qfvbq_1d(ji) = qfvbq_1d(ji) + & 714 715 709 fsup(ji) + ( 1.0 - zihgnew ) * & 710 focea(ji) * a_i_b(ji) * rdt_ice 716 711 717 712 zihic = 1.0 - MAX( zzero , SIGN( zone , -ht_i_b(ji) ) ) … … 719 714 fscbq_1d(ji) = a_i_b(ji) * fstbif_1d(ji) 720 715 qldif_1d(ji) = qldif_1d(ji) & 721 722 723 716 + fsup(ji) + ( 1.0 - zihgnew ) * focea(ji) * a_i_b(ji) & 717 * rdt_ice & 718 + ( 1.0 - zihic ) * fscbq_1d(ji) * rdt_ice 724 719 END DO ! ji 725 720 … … 743 738 ht_i_b(ji) = zhgnew(ji) 744 739 END DO ! ji 745 !746 !------------------------------------------------------------------------------|747 ! 6) Snow-Ice formation |748 !------------------------------------------------------------------------------|749 !740 ! 741 !------------------------------------------------------------------------------| 742 ! 6) Snow-Ice formation | 743 !------------------------------------------------------------------------------| 744 ! 750 745 ! When snow load excesses Archimede's limit, snow-ice interface goes down 751 746 ! under sea-level, flooding of seawater transforms snow into ice … … 754 749 755 750 dh_snowice(ji) = MAX(zzero,(rhosn*ht_s_b(ji)+(rhoic-rau0) & 756 751 * ht_i_b(ji))/(rhosn+rau0-rhoic)) 757 752 zhgnew(ji) = MAX(zhgnew(ji),zhgnew(ji)+dh_snowice(ji)) 758 753 zhnnew = MIN(ht_s_b(ji),ht_s_b(ji)-dh_snowice(ji)) 759 754 760 ! Changes in ice volume and ice mass.755 ! Changes in ice volume and ice mass. 761 756 dvnbq_1d(ji) = a_i_b(ji) * (zhgnew(ji)-ht_i_b(ji)) 762 757 dmgwi_1d(ji) = dmgwi_1d(ji) + a_i_b(ji) & 763 758 *(ht_s_b(ji)-zhnnew)*rhosn 764 759 765 760 rdmicif_1d(ji) = rdmicif_1d(ji) + a_i_b(ji) & 766 761 * ( zhgnew(ji) - ht_i_b(ji) )*rhoic 767 762 rdmsnif_1d(ji) = rdmsnif_1d(ji) + a_i_b(ji) & 768 769 770 ! Equivalent salt flux (1) Snow-ice formation component771 ! -----------------------------------------------------763 * ( zhnnew - ht_s_b(ji) )*rhosn 764 765 ! Equivalent salt flux (1) Snow-ice formation component 766 ! ----------------------------------------------------- 772 767 zji = MOD( npb(ji) - 1, jpi ) + 1 773 768 zjj = ( npb(ji) - 1 ) / jpi + 1 774 769 775 770 zsm_snowice = ( rhoic - rhosn ) / rhoic * & 776 771 sss_m(zji,zjj) 777 772 778 773 IF ( num_sal .NE. 2 ) zsm_snowice = sm_i_b(ji) 779 774 780 775 IF ( num_sal .NE. 4 ) & 781 fseqv_1d(ji) = fseqv_1d(ji) + &782 783 784 776 fseqv_1d(ji) = fseqv_1d(ji) + & 777 ( sss_m(zji,zjj) - zsm_snowice ) * & 778 a_i_b(ji) * & 779 ( zhgnew(ji) - ht_i_b(ji) ) * rhoic / rdt_ice 785 780 ! new lines 786 781 IF ( num_sal .EQ. 4 ) & 787 fseqv_1d(ji) = fseqv_1d(ji) + &788 789 790 782 fseqv_1d(ji) = fseqv_1d(ji) + & 783 ( sss_m(zji,zjj) - bulk_sal ) * & 784 a_i_b(ji) * & 785 ( zhgnew(ji) - ht_i_b(ji) ) * rhoic / rdt_ice 791 786 792 787 ! entrapment during snow ice formation 793 788 i_ice_switch = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_i_b(ji) + 1.0e-6 ) ) 794 789 isnowic = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - dh_snowice(ji) ) ) * & 795 790 i_ice_switch 796 791 IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) & 797 798 799 800 801 802 ! Actualize new snow and ice thickness.792 dsm_i_si_1d(ji) = ( zsm_snowice*dh_snowice(ji) & 793 + sm_i_b(ji) * ht_i_b(ji) & 794 / MAX( ht_i_b(ji) + dh_snowice(ji), zeps) & 795 - sm_i_b(ji) ) * isnowic 796 797 ! Actualize new snow and ice thickness. 803 798 ht_s_b(ji) = zhnnew 804 799 ht_i_b(ji) = zhgnew(ji) … … 811 806 zjj = ( npb(ji) - 1 ) / jpi + 1 812 807 diag_sni_gr(zji,zjj) = diag_sni_gr(zji,zjj) + dh_snowice(ji)*a_i_b(ji) / & 813 808 rdt_ice 814 809 815 810 END DO !ji 816 811 817 812 END SUBROUTINE lim_thd_dh 818 813 #else 819 814 !!====================================================================== … … 825 820 END SUBROUTINE lim_thd_dh 826 821 #endif 827 822 END MODULE limthd_dh
Note: See TracChangeset
for help on using the changeset viewer.
