Changeset 1572 for trunk/NEMO/LIM_SRC_3/limthd_dh.F90
- Timestamp:
- 2009-08-03T16:53:15+02:00 (15 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/LIM_SRC_3/limthd_dh.F90 (modified) (31 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/LIM_SRC_3/limthd_dh.F90
r1571 r1572 1 1 MODULE limthd_dh 2 !!====================================================================== 3 !! *** MODULE limthd_dh *** 4 !! LIM-3 : thermodynamic growth and decay of the ice 5 !!====================================================================== 6 !! History : LIM ! 2003-05 (M. Vancoppenolle) Original code in 1D 7 !! ! 2005-06 (M. Vancoppenolle) 3D version 8 !! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in rdmsnif & rdmicif 9 !!---------------------------------------------------------------------- 2 10 #if defined key_lim3 3 11 !!---------------------------------------------------------------------- 4 12 !! 'key_lim3' LIM3 sea-ice model 5 13 !!---------------------------------------------------------------------- 6 !!====================================================================== 7 !! *** MODULE limthd_dh *** 8 !! thermodynamic growth and decay of the ice 9 !!====================================================================== 10 14 !! lim_thd_dh : vertical accr./abl. and lateral ablation of sea ice 11 15 !!---------------------------------------------------------------------- 12 !! lim_thd_dh : vertical accr./abl. and lateral ablation of sea ice13 !! * Modules used14 15 16 USE par_oce ! ocean parameters 16 17 USE phycst ! physical constants (OCE directory) … … 27 28 PRIVATE 28 29 29 !! * Routine accessibility 30 PUBLIC lim_thd_dh ! called by lim_thd 31 32 !! * Module variables 33 REAL(wp) :: & ! constant values 34 epsi20 = 1e-20 , & 35 epsi13 = 1e-13 , & 36 epsi16 = 1e-16 , & 37 zzero = 0.e0 , & 38 zone = 1.e0 30 PUBLIC lim_thd_dh ! called by lim_thd 31 32 REAL(wp) :: epsi20 = 1e-20 ! constant values 33 REAL(wp) :: epsi13 = 1e-13 ! 34 REAL(wp) :: epsi16 = 1e-16 ! 35 REAL(wp) :: zzero = 0.e0 ! 36 REAL(wp) :: zone = 1.e0 ! 39 37 40 38 !!---------------------------------------------------------------------- 41 !! LIM 3.0, UCL-LOCEAN-IPSL (2005)39 !! NEMO/LIM 3.2, UCL-LOCEAN-IPSL (2009) 42 40 !! $Id$ 43 41 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 49 47 !!------------------------------------------------------------------ 50 48 !! *** ROUTINE lim_thd_dh *** 49 !! 50 !! ** Purpose : determines variations of ice and snow thicknesses. 51 !! 52 !! ** Method : Ice/Snow surface melting arises from imbalance in surface fluxes 53 !! Bottom accretion/ablation arises from flux budget 54 !! Snow thickness can increase by precipitation and decrease by sublimation 55 !! If snow load excesses Archmiede limit, snow-ice is formed by 56 !! the flooding of sea-water in the snow 57 !! 58 !! 1) Compute available flux of heat for surface ablation 59 !! 2) Compute snow and sea ice enthalpies 60 !! 3) Surface ablation and sublimation 61 !! 4) Bottom accretion/ablation 62 !! 5) Case of Total ablation 63 !! 6) Snow ice formation 64 !! 65 !! References : Bitz and Lipscomb, 1999, J. Geophys. Res. 66 !! Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646 67 !! Vancoppenolle, Fichefet and Bitz, 2005, Geophys. Res. Let. 68 !! Vancoppenolle et al.,2009, Ocean Modelling 51 69 !!------------------------------------------------------------------ 52 !! ** Purpose : 53 !! This routine determines variations of ice and snow thicknesses. 54 !! ** Method : 55 !! Ice/Snow surface melting arises from imbalance in surface fluxes 56 !! Bottom accretion/ablation arises from flux budget 57 !! Snow thickness can increase by precipitation and decrease by 58 !! sublimation 59 !! If snow load excesses Archmiede limit, snow-ice is formed by 60 !! the flooding of sea-water in the snow 61 !! ** Steps 62 !! 1) Compute available flux of heat for surface ablation 63 !! 2) Compute snow and sea ice enthalpies 64 !! 3) Surface ablation and sublimation 65 !! 4) Bottom accretion/ablation 66 !! 5) Case of Total ablation 67 !! 6) Snow ice formation 68 !! 69 !! ** Arguments 70 !! 71 !! ** Inputs / Outputs 72 !! 73 !! ** External 74 !! 75 !! ** References : Bitz and Lipscomb, JGR 99 76 !! Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646 77 !! Vancoppenolle, Fichefet and Bitz, GRL 2005 78 !! Vancoppenolle et al., OM08 79 !! 80 !! ** History : 81 !! original code 01-04 (LIM) 82 !! original routine 83 !! (05-2003) M. Vancoppenolle, Louvain-La-Neuve, Belgium 84 !! (06/07-2005) 3D version 85 !! (03-2008) Clean code 86 !! 70 INTEGER , INTENT(in) :: kideb, kiut ! Start/End point on which the the computation is applied 71 INTEGER , INTENT(in) :: jl ! Thickness cateogry number 72 !! 73 INTEGER :: ji , jk ! dummy loop indices 74 INTEGER :: zji, zjj ! 2D corresponding indices to ji 75 INTEGER :: isnow ! switch for presence (1) or absence (0) of snow 76 INTEGER :: isnowic ! snow ice formation not 77 INTEGER :: i_ice_switch ! ice thickness above a certain treshold or not 78 INTEGER :: iter 79 80 REAL(wp) :: zzfmass_i, zzfmass_s ! temporary scalar 81 REAL(wp) :: zhsnew, zihgnew, ztmelts ! temporary scalar 82 REAL(wp) :: zhn, zdhcf, zdhbf, zhni, zhnfi, zihg ! 83 REAL(wp) :: zdhnm, zhnnew, zhisn, zihic ! 84 REAL(wp) :: zfracs ! fractionation coefficient for bottom salt entrapment 85 REAL(wp) :: zds ! increment of bottom ice salinity 86 REAL(wp) :: zcoeff ! dummy argument for snowfall partitioning over ice and leads 87 REAL(wp) :: zsm_snowice ! snow-ice salinity 88 REAL(wp) :: zswi1 ! switch for computation of bottom salinity 89 REAL(wp) :: zswi12 ! switch for computation of bottom salinity 90 REAL(wp) :: zswi2 ! switch for computation of bottom salinity 91 REAL(wp) :: zgrr ! bottom growth rate 92 REAL(wp) :: ztform ! bottom formation temperature 93 94 REAL(wp), DIMENSION(jpij) :: zh_i ! ice layer thickness 95 REAL(wp), DIMENSION(jpij) :: zh_s ! snow layer thickness 96 REAL(wp), DIMENSION(jpij) :: ztfs ! melting point 97 REAL(wp), DIMENSION(jpij) :: zhsold ! old snow thickness 98 REAL(wp), DIMENSION(jpij) :: zqprec ! energy of fallen snow 99 REAL(wp), DIMENSION(jpij) :: zqfont_su ! incoming, remaining surface energy 100 REAL(wp), DIMENSION(jpij) :: zqfont_bo ! incoming, bottom energy 101 REAL(wp), DIMENSION(jpij) :: z_f_surf ! surface heat for ablation 102 REAL(wp), DIMENSION(jpij) :: zhgnew ! new ice thickness 103 REAL(wp), DIMENSION(jpij) :: zfmass_i ! 104 105 REAL(wp), DIMENSION(jpij) :: zdh_s_mel ! snow melt 106 REAL(wp), DIMENSION(jpij) :: zdh_s_pre ! snow precipitation 107 REAL(wp), DIMENSION(jpij) :: zdh_s_sub ! snow sublimation 108 REAL(wp), DIMENSION(jpij) :: zfsalt_melt ! salt flux due to ice melt 109 110 REAL(wp) , DIMENSION(jpij,jkmax) :: zdeltah 111 112 ! Pathological cases 113 REAL(wp), DIMENSION(jpij) :: zfdt_init ! total incoming heat for ice melt 114 REAL(wp), DIMENSION(jpij) :: zfdt_final ! total remaing heat for ice melt 115 REAL(wp), DIMENSION(jpij) :: zqt_i ! total ice heat content 116 REAL(wp), DIMENSION(jpij) :: zqt_s ! total snow heat content 117 REAL(wp), DIMENSION(jpij) :: zqt_dummy ! dummy heat content 118 119 REAL(wp), DIMENSION(jpij,jkmax) :: zqt_i_lay ! total ice heat content 120 121 ! Heat conservation 122 INTEGER :: num_iter_max, numce_dh 123 REAL(wp) :: meance_dh 124 INTEGER , DIMENSION(jpij) :: innermelt 125 REAL(wp), DIMENSION(jpij) :: zfbase, zdq_i 87 126 !!------------------------------------------------------------------ 88 !! * Arguments89 INTEGER , INTENT (IN) :: &90 kideb , & !: Start point on which the the computation is applied91 kiut , & !: End point on which the the computation is applied92 jl !: Thickness cateogry number93 94 !! * Local variables95 INTEGER :: &96 ji , & !: space index97 jk , & !: ice layer index98 isnow , & !: switch for presence (1) or absence (0) of snow99 zji , & !: 2D corresponding indices to ji100 zjj , &101 isnowic , & !: snow ice formation not102 i_ice_switch , & !: ice thickness above a certain treshold or not103 iter104 105 REAL(wp) :: &106 zzfmass_i , &107 zzfmass_s , &108 zhsnew , & !: new snow thickness109 zihgnew , & !: switch for total ablation110 ztmelts , & !: melting point111 zhn , &112 zdhcf , &113 zdhbf , &114 zhni , &115 zhnfi , &116 zihg , &117 zdhnm , &118 zhnnew , &119 zeps = 1.0e-13, &120 zhisn , &121 zfracs , & !: fractionation coefficient for bottom salt122 !: entrapment123 zds , & !: increment of bottom ice salinity124 zcoeff , & !: dummy argument for snowfall partitioning125 !: over ice and leads126 zsm_snowice, & !: snow-ice salinity127 zswi1 , & !: switch for computation of bottom salinity128 zswi12 , & !: switch for computation of bottom salinity129 zswi2 , & !: switch for computation of bottom salinity130 zgrr , & !: bottom growth rate131 zihic , & !:132 ztform !: bottom formation temperature133 134 REAL(wp) , DIMENSION(jpij) :: &135 zh_i , & ! ice layer thickness136 zh_s , & ! snow layer thickness137 ztfs , & ! melting point138 zhsold , & ! old snow thickness139 zqprec , & !: energy of fallen snow140 zqfont_su , & ! incoming, remaining surface energy141 zqfont_bo ! incoming, bottom energy142 143 REAL(wp) , DIMENSION(jpij) :: &144 z_f_surf, & ! surface heat for ablation145 zhgnew , & ! new ice thickness146 zfmass_i147 148 REAL(wp), DIMENSION(jpij) :: &149 zdh_s_mel , & ! snow melt150 zdh_s_pre , & ! snow precipitation151 zdh_s_sub , & ! snow sublimation152 zfsalt_melt ! salt flux due to ice melt153 154 REAL(wp) , DIMENSION(jpij,jkmax) :: &155 zdeltah156 157 ! Pathological cases158 REAL(wp), DIMENSION(jpij) :: &159 zfdt_init , & !: total incoming heat for ice melt160 zfdt_final , & !: total remaing heat for ice melt161 zqt_i , & !: total ice heat content162 zqt_s , & !: total snow heat content163 zqt_dummy !: dummy heat content164 165 REAL(wp), DIMENSION(jpij,jkmax) :: &166 zqt_i_lay !: total ice heat content167 168 ! Heat conservation169 REAL(wp), DIMENSION(jpij) :: &170 zfbase, &171 zdq_i172 173 INTEGER, DIMENSION(jpij) :: &174 innermelt175 176 REAL(wp) :: &177 meance_dh178 179 INTEGER :: &180 num_iter_max, &181 numce_dh182 127 183 128 zfsalt_melt(:) = 0.0 … … 198 143 isnow = INT( 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_s_b(ji) ) ) ) 199 144 ztfs(ji) = isnow * rtt + ( 1.0 - isnow ) * rtt 200 z_f_surf(ji) = qnsr_ice_1d(ji) + ( 1.0 - i0(ji) ) * & 201 qsr_ice_1d(ji) - fc_su(ji) 202 z_f_surf(ji) = MAX( zzero , z_f_surf(ji) ) * & 203 MAX( zzero , SIGN( zone , t_su_b(ji) - ztfs(ji) ) ) 204 zfdt_init(ji) = ( z_f_surf(ji) + & 205 MAX( fbif_1d(ji) + qlbbq_1d(ji) + fc_bo_i(ji),0.0 ) ) & 206 * rdt_ice 145 z_f_surf(ji) = qnsr_ice_1d(ji) + ( 1.0 - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 146 z_f_surf(ji) = MAX( zzero , z_f_surf(ji) ) * MAX( zzero , SIGN( zone , t_su_b(ji) - ztfs(ji) ) ) 147 zfdt_init(ji) = ( z_f_surf(ji) + MAX( fbif_1d(ji) + qlbbq_1d(ji) + fc_bo_i(ji),0.0 ) ) * rdt_ice 207 148 END DO ! ji 208 149 … … 226 167 DO jk = 1, nlay_s 227 168 DO ji = kideb,kiut 228 zqt_s(ji) 169 zqt_s(ji) = zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s 229 170 END DO 230 171 END DO … … 235 176 DO ji = kideb,kiut 236 177 zqt_i(ji) = zqt_i(ji) + q_i_b(ji,jk) * ht_i_b(ji) / nlay_i 237 zqt_i_lay(ji,jk) = q_i_b(ji,jk) * ht_i_b(ji) / nlay_i178 zqt_i_lay(ji,jk) = q_i_b(ji,jk) * ht_i_b(ji) / nlay_i 238 179 END DO 239 180 END DO … … 267 208 DO ji = kideb, kiut 268 209 ! tatm_ice is now in K 269 zqprec(ji) = rhosn * ( cpic * ( rtt - tatm_ice_1d(ji) ) + lfus ) 270 zqfont_su(ji) = z_f_surf(ji) * rdt_ice 271 zdeltah(ji,1) = MIN( 0.0 , - zqfont_su(ji) / MAX( zqprec(ji) , epsi13 ) ) 272 zqfont_su(ji) = MAX( 0.0 , - zdh_s_pre(ji) - zdeltah(ji,1) ) * & 273 zqprec(ji) 274 zdeltah(ji,1) = MAX( - zdh_s_pre(ji) , zdeltah(ji,1) ) 275 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,1) 210 zqprec (ji) = rhosn * ( cpic * ( rtt - tatm_ice_1d(ji) ) + lfus ) 211 zqfont_su(ji) = z_f_surf(ji) * rdt_ice 212 zdeltah (ji,1) = MIN( 0.e0 , - zqfont_su(ji) / MAX( zqprec(ji) , epsi13 ) ) 213 zqfont_su(ji) = MAX( 0.e0 , - zdh_s_pre(ji) - zdeltah(ji,1) ) * zqprec(ji) 214 zdeltah (ji,1) = MAX( - zdh_s_pre(ji) , zdeltah(ji,1) ) 215 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,1) 276 216 ! heat conservation 277 qt_s_in(ji,jl) = qt_s_in(ji,jl) + zqprec(ji) * zdh_s_pre(ji)278 zqt_s (ji) = zqt_s(ji)+ zqprec(ji) * zdh_s_pre(ji)279 zqt_s (ji) = MAX ( zqt_s(ji) - zqfont_su(ji) , 0.0 )217 qt_s_in(ji,jl) = qt_s_in(ji,jl) + zqprec(ji) * zdh_s_pre(ji) 218 zqt_s (ji) = zqt_s (ji) + zqprec(ji) * zdh_s_pre(ji) 219 zqt_s (ji) = MAX( zqt_s(ji) - zqfont_su(ji) , 0.e0 ) 280 220 END DO 281 221 … … 284 224 DO jk = 1, nlay_s 285 225 DO ji = kideb, kiut 286 zdeltah(ji,jk) = - zqfont_su(ji) / q_s_b(ji,jk) 287 zqfont_su(ji) = MAX( 0.0 , - zh_s(ji) - zdeltah(ji,jk) ) * & 288 q_s_b(ji,jk) 289 zdeltah(ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) ) 290 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) !resulting melt of snow 226 zdeltah (ji,jk) = - zqfont_su(ji) / q_s_b(ji,jk) 227 zqfont_su(ji) = MAX( 0.0 , - zh_s(ji) - zdeltah(ji,jk) ) * q_s_b(ji,jk) 228 zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) ) 229 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) ! resulting melt of snow 291 230 END DO 292 231 END DO … … 302 241 ht_s_b(ji) = MAX( zzero , zhsnew ) 303 242 ! Volume and mass variations of snow 304 dvsbq_1d(ji) = a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji) & 305 - zdh_s_mel(ji) ) 306 dvsbq_1d(ji) = MIN( zzero, dvsbq_1d(ji) ) 243 dvsbq_1d (ji) = a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji) - zdh_s_mel(ji) ) 244 dvsbq_1d (ji) = MIN( zzero, dvsbq_1d(ji) ) 307 245 rdmsnif_1d(ji) = rdmsnif_1d(ji) + rhosn * dvsbq_1d(ji) 308 246 END DO ! ji … … 312 250 !-------------------------- 313 251 DO ji = kideb, kiut 314 dh_i_surf(ji) = 0.0 315 ! For heat conservation test 316 z_f_surf(ji) = zqfont_su(ji) / rdt_ice ! heat conservation test 317 zdq_i(ji) = 0.0 252 dh_i_surf(ji) = 0.e0 253 z_f_surf (ji) = zqfont_su(ji) / rdt_ice ! heat conservation test 254 zdq_i (ji) = 0.e0 318 255 END DO ! ji 319 256 320 257 DO jk = 1, nlay_i 321 258 DO ji = kideb, kiut 322 ! melt of layer jk 323 zdeltah(ji,jk) = - zqfont_su(ji) / q_i_b(ji,jk) 324 ! recompute heat available 325 zqfont_su(ji) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * & 326 q_i_b(ji,jk) 327 ! melt of layer jk cannot be higher than its thickness 328 zdeltah(ji,jk) = MAX( zdeltah(ji,jk) , - zh_i(ji) ) 329 ! update surface melt 330 dh_i_surf(ji) = dh_i_surf(ji) + zdeltah(ji,jk) 331 ! for energy conservation 332 zdq_i(ji) = zdq_i(ji) + zdeltah(ji,jk) * & 333 q_i_b(ji,jk) / rdt_ice 259 ! ! melt of layer jk 260 zdeltah (ji,jk) = - zqfont_su(ji) / q_i_b(ji,jk) 261 ! ! recompute heat available 262 zqfont_su(ji) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * q_i_b(ji,jk) 263 ! ! melt of layer jk cannot be higher than its thickness 264 zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_i(ji) ) 265 ! ! update surface melt 266 dh_i_surf(ji) = dh_i_surf(ji) + zdeltah(ji,jk) 267 ! ! for energy conservation 268 zdq_i (ji) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) / rdt_ice 269 ! 334 270 ! contribution to ice-ocean salt flux 335 zji = MOD( npb(ji) - 1, jpi ) + 1 336 zjj = ( npb(ji) - 1 ) / jpi + 1 337 zfsalt_melt(ji) = zfsalt_melt(ji) + & 338 ( sss_m(zji,zjj) - sm_i_b(ji) ) * & 339 a_i_b(ji) * & 340 MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice 341 END DO ! ji 342 END DO ! jk 343 344 !------------------- 345 ! Conservation test 346 !------------------- 347 IF ( con_i ) THEN 348 numce_dh = 0 349 meance_dh = 0.0 271 zji = MOD( npb(ji) - 1, jpi ) + 1 272 zjj = ( npb(ji) - 1 ) / jpi + 1 273 zfsalt_melt(ji) = zfsalt_melt(ji) + ( sss_m(zji,zjj) - sm_i_b(ji) ) * a_i_b(ji) & 274 & * MIN( zdeltah(ji,jk) , 0.e0 ) * rhoic / rdt_ice 275 END DO 276 END DO 277 278 ! !------------------- 279 IF( con_i ) THEN ! Conservation test 280 ! !------------------- 281 numce_dh = 0 282 meance_dh = 0.e0 350 283 DO ji = kideb, kiut 351 352 284 IF ( ( z_f_surf(ji) + zdq_i(ji) ) .GE. 1.0e-3 ) THEN 353 285 numce_dh = numce_dh + 1 354 286 meance_dh = meance_dh + z_f_surf(ji) + zdq_i(ji) 355 287 ENDIF 356 357 IF ( z_f_surf(ji) + zdq_i(ji) .GE. 1.0e-3 ) THEN! 288 IF( z_f_surf(ji) + zdq_i(ji) .GE. 1.0e-3 ) THEN! 358 289 WRITE(numout,*) ' ALERTE heat loss for surface melt ' 359 290 WRITE(numout,*) ' zji, zjj, jl :', zji, zjj, jl … … 368 299 WRITE(numout,*) ' sss_m : ', sss_m(zji,zjj) 369 300 ENDIF 370 371 END DO ! ji 372 373 IF ( numce_dh .GT. 0 ) meance_dh = meance_dh / numce_dh 301 END DO 302 ! 303 IF( numce_dh > 0 ) meance_dh = meance_dh / numce_dh 374 304 WRITE(numout,*) ' Error report - Category : ', jl 375 305 WRITE(numout,*) ' ~~~~~~~~~~~~ ' 376 306 WRITE(numout,*) ' Number of points where there is sur. me. error : ', numce_dh 377 307 WRITE(numout,*) ' Mean basal growth error on error points : ', meance_dh 378 379 ENDIF ! con_i308 ! 309 ENDIF 380 310 381 311 !---------------------- … … 386 316 ! if qla is positive (upwards), heat goes to the atmosphere, therefore 387 317 ! snow sublimates, if qla is negative (downwards), snow condensates 388 zdh_s_sub(ji) = - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice389 dh_s_tot (ji) = dh_s_tot(ji) + zdh_s_sub(ji)390 zdhcf = ht_s_b(ji) + zdh_s_sub(ji)391 ht_s_b (ji)= MAX( zzero , zdhcf )318 zdh_s_sub(ji) = - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice 319 dh_s_tot (ji) = dh_s_tot(ji) + zdh_s_sub(ji) 320 zdhcf = ht_s_b(ji) + zdh_s_sub(ji) 321 ht_s_b (ji) = MAX( zzero , zdhcf ) 392 322 ! we recompute dh_s_tot 393 dh_s_tot (ji) = ht_s_b(ji) - zhsold(ji)394 qt_s_in (ji,jl) = qt_s_in(ji,jl) + zdh_s_sub(ji)*q_s_b(ji,1)395 END DO !ji396 397 zqt_dummy(:) = 0. 0323 dh_s_tot (ji) = ht_s_b(ji) - zhsold(ji) 324 qt_s_in (ji,jl) = qt_s_in(ji,jl) + zdh_s_sub(ji)*q_s_b(ji,1) 325 END DO 326 327 zqt_dummy(:) = 0.e0 398 328 DO jk = 1, nlay_s 399 329 DO ji = kideb,kiut 400 q_s_b(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) 401 ! heat conservation 402 zqt_dummy(ji) = zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s 403 END DO 404 END DO 405 406 DO jk = 1, nlay_s !n 407 DO ji = kideb, kiut !n 408 ! In case of disparition of the snow, we have to update the snow 409 ! temperatures 330 q_s_b (ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) 331 zqt_dummy(ji) = zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s ! heat conservation 332 END DO 333 END DO 334 335 DO jk = 1, nlay_s 336 DO ji = kideb, kiut 337 ! In case of disparition of the snow, we have to update the snow temperatures 410 338 zhisn = MAX( zzero , SIGN( zone, - ht_s_b(ji) ) ) 411 339 t_s_b(ji,jk) = ( 1.0 - zhisn ) * t_s_b(ji,jk) + zhisn * rtt … … 428 356 ! 4.1 Basal growth - (a) salinity not varying in time 429 357 !----------------------------------------------------- 430 IF ( ( num_sal .NE. 2 ) .AND. ( num_sal .NE. 4 )) THEN358 IF( num_sal /= 2 .AND. num_sal /= 4 ) THEN 431 359 DO ji = kideb, kiut 432 IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT.0.0 ) THEN360 IF( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) < 0.0 ) THEN 433 361 s_i_new(ji) = sm_i_b(ji) 434 362 ! Melting point in K … … 437 365 ztform = t_i_b(ji,nlay_i) ! t_bo_b crashes in the 438 366 ! Baltic 439 q_i_b(ji,nlay_i+1) = rhoic * & 440 ( cpic * ( ztmelts - ztform ) & 441 + lfus * ( 1.0 - ( ztmelts - rtt ) / & 442 ( ztform - rtt ) ) & 443 - rcp * ( ztmelts-rtt ) ) 367 q_i_b(ji,nlay_i+1) = rhoic * ( cpic * ( ztmelts - ztform ) & 368 & + lfus * ( 1.0 - ( ztmelts - rtt ) / ( ztform - rtt ) ) & 369 & - rcp * ( ztmelts - rtt ) ) 444 370 ! Basal growth rate = - F*dt / q 445 dh_i_bott(ji) = - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + & 446 qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 447 ENDIF ! heat budget 448 END DO ! ji 449 ENDIF ! num_sal 371 dh_i_bott(ji) = - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 372 ENDIF 373 END DO 374 ENDIF 450 375 451 376 !------------------------------------------------- 452 377 ! 4.1 Basal growth - (b) salinity varying in time 453 378 !------------------------------------------------- 454 IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 )) THEN379 IF( num_sal == 2 .OR. num_sal == 4 ) THEN 455 380 ! the growth rate (dh_i_bott) is function of the new ice 456 381 ! heat content (q_i_b(nlay_i+1)). q_i_b depends on the new ice … … 461 386 ! Initial value (tested 1D, can be anything between 1 and 20) 462 387 num_iter_max = 4 463 s_i_new(:) = 4.0388 s_i_new(:) = 4.0 464 389 465 390 ! Iterative procedure 466 391 DO iter = 1, num_iter_max 467 392 DO ji = kideb, kiut 468 IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0 ) THEN393 IF( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) < 0.e0 ) THEN 469 394 zji = MOD( npb(ji) - 1, jpi ) + 1 470 395 zjj = ( npb(ji) - 1 ) / jpi + 1 … … 472 397 ztmelts = - tmut * s_i_new(ji) + rtt 473 398 ! New ice heat content (Bitz and Lipscomb, 1999) 474 q_i_b(ji,nlay_i+1) = rhoic * & 475 ( cpic * ( ztmelts - t_bo_b(ji) ) & 476 + lfus * ( 1.0 - ( ztmelts - rtt ) / & 477 ( t_bo_b(ji) - rtt ) ) & 478 - rcp * ( ztmelts-rtt ) ) 399 q_i_b(ji,nlay_i+1) = rhoic * ( cpic * ( ztmelts - t_bo_b(ji) ) & 400 & + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) ) & 401 & - rcp * ( ztmelts-rtt ) ) 479 402 ! Bottom growth rate = - F*dt / q 480 dh_i_bott(ji) = - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) & 481 + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 403 dh_i_bott(ji) = - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 482 404 ! New ice salinity ( Cox and Weeks, JGR, 1988 ) 483 405 ! zswi2 (1) if dh_i_bott/rdt .GT. 3.6e-7 484 406 ! zswi12 (1) if dh_i_bott/rdt .LT. 3.6e-7 and .GT. 2.0e-8 485 407 ! zswi1 (1) if dh_i_bott/rdt .LT. 2.0e-8 486 zgrr = MIN( 1.0e-3, MAX ( dh_i_bott(ji) / rdt_ice , zeps) )408 zgrr = MIN( 1.0e-3, MAX ( dh_i_bott(ji) / rdt_ice , epsi13 ) ) 487 409 zswi2 = MAX( zzero , SIGN( zone , zgrr - 3.6e-7 ) ) 488 410 zswi12 = MAX( zzero , SIGN( zone , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 ) 489 411 zswi1 = 1. - zswi2 * zswi12 490 zfracs = zswi1 * 0.12 + & 491 zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) ) + & 492 zswi2 * 0.26 / & 493 ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) 494 zds = zfracs*sss_m(zji,zjj) - s_i_new(ji) 412 zfracs = zswi1 * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) ) & 413 & + zswi2 * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) 414 zds = zfracs * sss_m(zji,zjj) - s_i_new(ji) 495 415 s_i_new(ji) = zfracs * sss_m(zji,zjj) 496 416 ENDIF ! fc_bo_i … … 500 420 ! Final values 501 421 DO ji = kideb, kiut 502 IF 422 IF( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0 ) THEN 503 423 ! New ice salinity must not exceed 15 psu 504 424 s_i_new(ji) = MIN( s_i_new(ji), s_i_max ) … … 506 426 ztmelts = - tmut * s_i_new(ji) + rtt 507 427 ! New ice heat content (Bitz and Lipscomb, 1999) 508 q_i_b(ji,nlay_i+1) = rhoic * & 509 ( cpic * ( ztmelts - t_bo_b(ji) ) & 510 + lfus * ( 1.0 - ( ztmelts - rtt ) / & 511 ( t_bo_b(ji) - rtt ) ) & 512 - rcp * ( ztmelts-rtt ) ) 428 q_i_b(ji,nlay_i+1) = rhoic * ( cpic * ( ztmelts - t_bo_b(ji) ) & 429 & + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) ) & 430 & - rcp * ( ztmelts - rtt ) ) 513 431 ! Basal growth rate = - F*dt / q 514 dh_i_bott(ji) = - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + & 515 qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 432 dh_i_bott(ji) = - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 516 433 ! Salinity update 517 434 ! entrapment during bottom growth 518 dsm_i_se_1d(ji) = ( s_i_new(ji)*dh_i_bott(ji) + & 519 sm_i_b(ji)*ht_i_b(ji) ) / & 520 MAX( ht_i_b(ji) + dh_i_bott(ji) ,zeps ) & 521 - sm_i_b(ji) 435 dsm_i_se_1d(ji) = ( s_i_new(ji) * dh_i_bott(ji) + sm_i_b(ji) * ht_i_b(ji) ) & 436 & / MAX( ht_i_b(ji) + dh_i_bott(ji) ,epsi13 ) - sm_i_b(ji) 522 437 ENDIF ! heat budget 523 END DO ! ji524 ENDIF ! num_sal438 END DO 439 ENDIF 525 440 526 441 !---------------- … … 528 443 !---------------- 529 444 meance_dh = 0.0 530 numce_dh = 0445 numce_dh = 0 531 446 innermelt(:) = 0 532 447 533 448 DO ji = kideb, kiut 534 449 ! heat convergence at the surface > 0 535 IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0 ) THEN450 IF( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) >= 0.e0 ) THEN 536 451 537 452 s_i_new(ji) = s_i_b(ji,nlay_i) … … 539 454 540 455 zfbase(ji) = zqfont_bo(ji) / rdt_ice ! heat conservation test 541 zdq_i(ji) = 0. 0542 543 dh_i_bott(ji) = 0. 0456 zdq_i(ji) = 0.e0 457 458 dh_i_bott(ji) = 0.e0 544 459 ENDIF 545 460 END DO … … 555 470 ELSE ! normal ablation 556 471 zdeltah(ji,jk) = - zqfont_bo(ji) / q_i_b(ji,jk) 557 zqfont_bo(ji) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * & 558 q_i_b(ji,jk) 472 zqfont_bo(ji) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * q_i_b(ji,jk) 559 473 zdeltah(ji,jk) = MAX(zdeltah(ji,jk), - zh_i(ji) ) 560 474 dh_i_bott(ji) = dh_i_bott(ji) + zdeltah(ji,jk) 561 zdq_i(ji) = zdq_i(ji) + zdeltah(ji,jk) * & 562 q_i_b(ji,jk) / rdt_ice 475 zdq_i(ji) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) / rdt_ice 563 476 ! contribution to salt flux 564 477 zji = MOD( npb(ji) - 1, jpi ) + 1 565 478 zjj = ( npb(ji) - 1 ) / jpi + 1 566 zfsalt_melt(ji) = zfsalt_melt(ji) + & 567 ( sss_m(zji,zjj) - sm_i_b(ji) ) * & 568 a_i_b(ji) * & 569 MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice 479 zfsalt_melt(ji) = zfsalt_melt(ji) + ( sss_m(zji,zjj) - sm_i_b(ji) ) * a_i_b(ji) & 480 & * MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice 570 481 ENDIF 571 482 ENDIF … … 573 484 END DO ! jk 574 485 575 !------------------- 576 ! Conservation test 577 !------------------- 578 IF ( con_i ) THEN 486 ! !------------------- 487 IF( con_i ) THEN ! Conservation test 488 ! !------------------- 579 489 DO ji = kideb, kiut 580 IF ( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0 ) THEN581 IF ( ( zfbase(ji) + zdq_i(ji) ) .GE. 1.0e-3 ) THEN582 numce_dh = numce_dh + 1490 IF( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) >= 0.e0 ) THEN 491 IF( ( zfbase(ji) + zdq_i(ji) ) >= 1.e-3 ) THEN 492 numce_dh = numce_dh + 1 583 493 meance_dh = meance_dh + zfbase(ji) + zdq_i(ji) 584 494 ENDIF … … 598 508 WRITE(numout,*) ' innermelt : ', innermelt(ji) 599 509 ENDIF 600 ENDIF ! heat convergence at the surface 601 END DO ! ji 602 603 IF ( numce_dh .GT. 0 ) meance_dh = meance_dh / numce_dh 510 ENDIF 511 END DO 512 IF( numce_dh > 0 ) meance_dh = meance_dh / numce_dh 604 513 WRITE(numout,*) ' Number of points where there is bas. me. error : ', numce_dh 605 514 WRITE(numout,*) ' Mean basal melt error on error points : ', meance_dh 606 515 WRITE(numout,*) ' Remaining bottom heat : ', zqfont_bo(jiindex_1d) 607 608 ENDIF ! con_i516 ! 517 ENDIF 609 518 610 519 ! … … 618 527 619 528 DO ji = kideb, kiut 620 ! in a 1-category sea ice model, bottom ablation must not exceed hmelt (-0.15)621 zdhbf = dh_i_bott(ji)622 IF (jpl.EQ.1) zdhbf = MAX( hmelt , dh_i_bott(ji) )623 ! excessive energy is sent to lateral ablation624 fsup(ji) = rhoic*lfus * at_i_b(ji) / MAX( ( 1.0 - at_i_b(ji) ),epsi13) &625 * ( zdhbf - dh_i_bott(ji) ) / rdt_ice626 529 ! ! in a 1-category sea ice model, bottom ablation must not exceed hmelt (-0.15) 530 IF( jpl == 1 ) THEN ; zdhbf = MAX( hmelt , dh_i_bott(ji) ) 531 ELSE ; zdhbf = dh_i_bott(ji) 532 ENDIF 533 ! ! excessive energy is sent to lateral ablation 534 fsup (ji) = rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi13 ) & 535 & * ( zdhbf - dh_i_bott(ji) ) / rdt_ice 627 536 dh_i_bott(ji) = zdhbf 628 !since ice volume is only used for outputs, we keep it global for all categories 629 dvbbq_1d(ji) = a_i_b(ji)*dh_i_bott(ji) 630 !new ice thickness 631 zhgnew(ji) = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji) 632 633 ! diagnostic ( bottom ice growth ) 634 zji = MOD( npb(ji) - 1, jpi ) + 1 635 zjj = ( npb(ji) - 1 ) / jpi + 1 636 diag_bot_gr(zji,zjj) = diag_bot_gr(zji,zjj) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) & 637 / rdt_ice 638 diag_sur_me(zji,zjj) = diag_sur_me(zji,zjj) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) & 639 / rdt_ice 640 diag_bot_me(zji,zjj) = diag_bot_me(zji,zjj) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) & 641 / rdt_ice 537 ! !since ice volume is only used for outputs, we keep it global for all categories 538 dvbbq_1d (ji) = a_i_b(ji) * dh_i_bott(ji) 539 ! !new ice thickness 540 zhgnew (ji) = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji) 541 ! ! diagnostic ( bottom ice growth ) 542 zji = MOD( npb(ji) - 1, jpi ) + 1 543 zjj = ( npb(ji) - 1 ) / jpi + 1 544 diag_bot_gr(zji,zjj) = diag_bot_gr(zji,zjj) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) / rdt_ice 545 diag_sur_me(zji,zjj) = diag_sur_me(zji,zjj) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) / rdt_ice 546 diag_bot_me(zji,zjj) = diag_bot_me(zji,zjj) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) / rdt_ice 642 547 END DO 643 548 … … 653 558 zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice 654 559 ! 0 if no more ice 655 zhgnew (ji) = zihgnew * zhgnew(ji)! ice thickness is put to 0560 zhgnew (ji) = zihgnew * zhgnew(ji) ! ice thickness is put to 0 656 561 ! remaining heat 657 562 zfdt_final(ji) = ( 1.0 - zihgnew ) * ( zqfont_su(ji) + zqfont_bo(ji) ) 658 563 659 564 ! If snow remains, energy is used to melt snow 660 zhni = ht_s_b(ji)! snow depth at previous time step661 zihg 565 zhni = ht_s_b(ji) ! snow depth at previous time step 566 zihg = MAX( zzero , SIGN ( zone , - ht_s_b(ji) ) ) ! 0 if snow 662 567 663 568 ! energy of melting of remaining snow 664 zqt_s(ji) = ( 1. - zihg) * zqt_s(ji) / MAX( zhni, zeps ) 665 zdhnm = - ( 1. - zihg ) * ( 1. - zihgnew ) * ( zfdt_final(ji) / & 666 MAX( zqt_s(ji) , zeps ) ) 569 zqt_s(ji) = ( 1. - zihg ) * zqt_s(ji) / MAX( zhni, epsi13 ) 570 zdhnm = - ( 1. - zihg ) * ( 1. - zihgnew ) * zfdt_final(ji) / MAX( zqt_s(ji) , epsi13 ) 667 571 zhnfi = zhni + zdhnm 668 zfdt_final(ji) = MAX 572 zfdt_final(ji) = MAX( zfdt_final(ji) + zqt_s(ji) * zdhnm , 0.0 ) 669 573 ht_s_b(ji) = MAX( zzero , zhnfi ) 670 574 zqt_s(ji) = zqt_s(ji) * ht_s_b(ji) … … 672 576 ! Mass variations of ice and snow 673 577 !--------------------------------- 674 ! 578 ! ! mass variation of the jl category 675 579 zzfmass_s = - a_i_b(ji) * ( zhni - ht_s_b(ji) ) * rhosn ! snow 676 580 zzfmass_i = a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic ! ice … … 684 588 ! Remaining heat to the ocean 685 589 !--------------------------------- 686 ! focea is in W.m-2 * dt 687 focea(ji) = - zfdt_final(ji) / rdt_ice 590 focea(ji) = - zfdt_final(ji) / rdt_ice ! focea is in W.m-2 * dt 688 591 689 592 END DO … … 695 598 !--------------------------- 696 599 DO ji = kideb, kiut 697 zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice600 zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice 698 601 699 602 ! Salt flux 700 zji = MOD( npb(ji) - 1, jpi ) + 1 701 zjj = ( npb(ji) - 1 ) / jpi + 1 702 IF ( num_sal .NE. 4 ) & 703 fseqv_1d(ji) = fseqv_1d(ji) + zihgnew * zfsalt_melt(ji) + & 704 (1.0 - zihgnew) * zfmass_i(ji) * & 705 ( sss_m(zji,zjj) - sm_i_b(ji) ) / rdt_ice 603 zji = MOD( npb(ji) - 1, jpi ) + 1 604 zjj = ( npb(ji) - 1 ) / jpi + 1 706 605 ! new lines 707 IF ( num_sal .EQ. 4 ) & 708 fseqv_1d(ji) = fseqv_1d(ji) + zihgnew * zfsalt_melt(ji) + & 709 (1.0 - zihgnew) * zfmass_i(ji) * & 710 ( sss_m(zji,zjj) - bulk_sal ) / rdt_ice 606 IF( num_sal == 4 ) THEN 607 fseqv_1d(ji) = fseqv_1d(ji) + zihgnew * zfsalt_melt(ji) & 608 & + (1.0 - zihgnew) * zfmass_i(ji) * ( sss_m(zji,zjj) - bulk_sal ) / rdt_ice 609 ELSE 610 fseqv_1d(ji) = fseqv_1d(ji) + zihgnew * zfsalt_melt(ji) & 611 & + (1.0 - zihgnew) * zfmass_i(ji) * ( sss_m(zji,zjj) - sm_i_b(ji) ) / rdt_ice 612 ENDIF 711 613 ! Heat flux 712 614 ! excessive bottom ablation energy (fsup) - 0 except if jpl = 1 713 615 ! excessive total ablation energy (focea) sent to the ocean 714 qfvbq_1d(ji) = qfvbq_1d(ji) + & 715 fsup(ji) + ( 1.0 - zihgnew ) * & 716 focea(ji) * a_i_b(ji) * rdt_ice 616 qfvbq_1d(ji) = qfvbq_1d(ji) + fsup(ji) + ( 1.0 - zihgnew ) * focea(ji) * a_i_b(ji) * rdt_ice 717 617 718 618 zihic = 1.0 - MAX( zzero , SIGN( zone , -ht_i_b(ji) ) ) 719 619 ! equals 0 if ht_i = 0, 1 if ht_i gt 0 720 620 fscbq_1d(ji) = a_i_b(ji) * fstbif_1d(ji) 721 qldif_1d(ji) = qldif_1d(ji) & 722 + fsup(ji) + ( 1.0 - zihgnew ) * focea(ji) * a_i_b(ji) & 723 * rdt_ice & 724 + ( 1.0 - zihic ) * fscbq_1d(ji) * rdt_ice 621 qldif_1d(ji) = qldif_1d(ji) + fsup(ji) + ( 1.0 - zihgnew ) * focea(ji) * a_i_b(ji) * rdt_ice & 622 & + ( 1.0 - zihic ) * fscbq_1d(ji) * rdt_ice 725 623 END DO ! ji 726 624 … … 729 627 !------------------------------------------- 730 628 DO ji = kideb, kiut 731 zihgnew 732 t_su_b(ji) 629 zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) 630 t_su_b(ji) = zihgnew * t_su_b(ji) + ( 1.0 - zihgnew ) * rtt 733 631 END DO ! ji 734 632 735 633 DO jk = 1, nlay_i 736 634 DO ji = kideb, kiut 737 zihgnew 738 t_i_b(ji,jk) 739 q_i_b(ji,jk) 635 zihgnew = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) 636 t_i_b(ji,jk) = zihgnew * t_i_b(ji,jk) + ( 1.0 - zihgnew ) * rtt 637 q_i_b(ji,jk) = zihgnew * q_i_b(ji,jk) 740 638 END DO 741 639 END DO ! ji … … 748 646 ! 6) Snow-Ice formation | 749 647 !------------------------------------------------------------------------------| 750 ! 751 ! When snow load excesses Archimede's limit, snow-ice interface goes down 752 ! under sea-level, flooding of seawater transforms snow into ice 753 ! dh_snowice is positive for the ice 754 DO ji = kideb, kiut 755 756 dh_snowice(ji) = MAX(zzero,(rhosn*ht_s_b(ji)+(rhoic-rau0) & 757 * ht_i_b(ji))/(rhosn+rau0-rhoic)) 758 zhgnew(ji) = MAX(zhgnew(ji),zhgnew(ji)+dh_snowice(ji)) 759 zhnnew = MIN(ht_s_b(ji),ht_s_b(ji)-dh_snowice(ji)) 648 ! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level, 649 ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice 650 DO ji = kideb, kiut 651 ! 652 dh_snowice(ji) = MAX( zzero , ( rhosn * ht_s_b(ji) + (rhoic-rau0) * ht_i_b(ji) ) / ( rhosn+rau0-rhoic ) ) 653 zhgnew(ji) = MAX( zhgnew(ji) , zhgnew(ji) + dh_snowice(ji) ) 654 zhnnew = MIN( ht_s_b(ji) , ht_s_b(ji) - dh_snowice(ji) ) 760 655 761 656 ! Changes in ice volume and ice mass. 762 dvnbq_1d(ji) = a_i_b(ji) * (zhgnew(ji)-ht_i_b(ji)) 763 dmgwi_1d(ji) = dmgwi_1d(ji) + a_i_b(ji) & 764 *(ht_s_b(ji)-zhnnew)*rhosn 765 766 rdmicif_1d(ji) = rdmicif_1d(ji) + a_i_b(ji) & 767 * ( zhgnew(ji) - ht_i_b(ji) )*rhoic 768 rdmsnif_1d(ji) = rdmsnif_1d(ji) + a_i_b(ji) & 769 * ( zhnnew - ht_s_b(ji) )*rhosn 657 dvnbq_1d (ji) = a_i_b(ji) * ( zhgnew(ji)-ht_i_b(ji) ) 658 dmgwi_1d (ji) = dmgwi_1d(ji) + a_i_b(ji) * ( ht_s_b(ji) - zhnnew ) * rhosn 659 660 rdmicif_1d(ji) = rdmicif_1d(ji) + a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic 661 rdmsnif_1d(ji) = rdmsnif_1d(ji) + a_i_b(ji) * ( zhnnew - ht_s_b(ji) ) * rhosn 770 662 771 663 ! Equivalent salt flux (1) Snow-ice formation component 772 664 ! ----------------------------------------------------- 773 zji = MOD( npb(ji) - 1, jpi ) + 1 774 zjj = ( npb(ji) - 1 ) / jpi + 1 775 776 zsm_snowice = ( rhoic - rhosn ) / rhoic * & 777 sss_m(zji,zjj) 778 779 IF ( num_sal .NE. 2 ) zsm_snowice = sm_i_b(ji) 780 781 IF ( num_sal .NE. 4 ) & 782 fseqv_1d(ji) = fseqv_1d(ji) + & 783 ( sss_m(zji,zjj) - zsm_snowice ) * & 784 a_i_b(ji) * & 785 ( zhgnew(ji) - ht_i_b(ji) ) * rhoic / rdt_ice 786 ! new lines 787 IF ( num_sal .EQ. 4 ) & 788 fseqv_1d(ji) = fseqv_1d(ji) + & 789 ( sss_m(zji,zjj) - bulk_sal ) * & 790 a_i_b(ji) * & 791 ( zhgnew(ji) - ht_i_b(ji) ) * rhoic / rdt_ice 792 665 zji = MOD( npb(ji) - 1, jpi ) + 1 666 zjj = ( npb(ji) - 1 ) / jpi + 1 667 668 IF( num_sal /= 2 ) THEN ; zsm_snowice = sm_i_b(ji) 669 ELSE ; zsm_snowice = ( rhoic - rhosn ) / rhoic * sss_m(zji,zjj) 670 ENDIF 671 IF( num_sal == 4 ) THEN 672 fseqv_1d(ji) = fseqv_1d(ji) + ( sss_m(zji,zjj) - bulk_sal ) * a_i_b(ji) & 673 & * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic / rdt_ice 674 ELSE 675 fseqv_1d(ji) = fseqv_1d(ji) + ( sss_m(zji,zjj) - zsm_snowice ) * a_i_b(ji) & 676 & * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic / rdt_ice 677 ENDIF 793 678 ! entrapment during snow ice formation 794 i_ice_switch = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_i_b(ji) + 1.0e-6 ) ) 795 isnowic = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - dh_snowice(ji) ) ) * & 796 i_ice_switch 797 IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) & 798 dsm_i_si_1d(ji) = ( zsm_snowice*dh_snowice(ji) & 799 + sm_i_b(ji) * ht_i_b(ji) & 800 / MAX( ht_i_b(ji) + dh_snowice(ji), zeps) & 801 - sm_i_b(ji) ) * isnowic 679 i_ice_switch = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - ht_i_b(ji) + 1.0e-6 ) ) 680 isnowic = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - dh_snowice(ji) ) ) * i_ice_switch 681 IF( num_sal == 2 .OR. num_sal == 4 ) & 682 dsm_i_si_1d(ji) = ( zsm_snowice*dh_snowice(ji) & 683 & + sm_i_b(ji) * ht_i_b(ji) / MAX( ht_i_b(ji) + dh_snowice(ji), epsi13) & 684 & - sm_i_b(ji) ) * isnowic 802 685 803 686 ! Actualize new snow and ice thickness. … … 806 689 807 690 ! Total ablation ! new lines added to debug 808 IF( ht_i_b(ji) .LE.0.0 )a_i_b(ji) = 0.0691 IF( ht_i_b(ji) <= 0.e0 ) a_i_b(ji) = 0.0 809 692 810 693 ! diagnostic ( snow ice growth ) 811 zji = MOD( npb(ji) - 1, jpi ) + 1 812 zjj = ( npb(ji) - 1 ) / jpi + 1 813 diag_sni_gr(zji,zjj) = diag_sni_gr(zji,zjj) + dh_snowice(ji)*a_i_b(ji) / & 814 rdt_ice 815 694 zji = MOD( npb(ji) - 1, jpi ) + 1 695 zjj = ( npb(ji) - 1 ) / jpi + 1 696 diag_sni_gr(zji,zjj) = diag_sni_gr(zji,zjj) + dh_snowice(ji)*a_i_b(ji) / rdt_ice 697 ! 816 698 END DO !ji 817 699 818 700 END SUBROUTINE lim_thd_dh 701 819 702 #else 820 !!====================================================================== 821 !! *** MODULE limthd_dh *** 822 !! no sea ice model 823 !!====================================================================== 703 !!---------------------------------------------------------------------- 704 !! Default option NO LIM3 sea-ice model 705 !!---------------------------------------------------------------------- 824 706 CONTAINS 825 707 SUBROUTINE lim_thd_dh ! Empty routine 826 708 END SUBROUTINE lim_thd_dh 827 709 #endif 710 711 !!====================================================================== 828 712 END MODULE limthd_dh
Note: See TracChangeset
for help on using the changeset viewer.
