- Timestamp:
- 2015-07-21T13:25:36+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r5134 r5621 29 29 PRIVATE 30 30 31 PUBLIC lim_thd_dh ! called by lim_thd 31 PUBLIC lim_thd_dh ! called by lim_thd 32 PUBLIC lim_thd_snwblow ! called in sbcblk/sbcclio/sbccpl and here 33 34 INTERFACE lim_thd_snwblow 35 MODULE PROCEDURE lim_thd_snwblow_1d, lim_thd_snwblow_2d 36 END INTERFACE 32 37 33 38 !!---------------------------------------------------------------------- … … 71 76 REAL(wp) :: zfdum 72 77 REAL(wp) :: zfracs ! fractionation coefficient for bottom salt entrapment 73 REAL(wp) :: zcoeff ! dummy argument for snowfall partitioning over ice and leads 74 REAL(wp) :: zs_snic ! snow-ice salinity 78 REAL(wp) :: zs_snic ! snow-ice salinity 75 79 REAL(wp) :: zswi1 ! switch for computation of bottom salinity 76 80 REAL(wp) :: zswi12 ! switch for computation of bottom salinity … … 86 90 REAL(wp) :: zsstK ! SST in Kelvin 87 91 88 REAL(wp), POINTER, DIMENSION(:) :: zh_s ! snow layer thickness89 92 REAL(wp), POINTER, DIMENSION(:) :: zqprec ! energy of fallen snow (J.m-3) 90 93 REAL(wp), POINTER, DIMENSION(:) :: zq_su ! heat for surface ablation (J.m-2) … … 92 95 REAL(wp), POINTER, DIMENSION(:) :: zq_rema ! remaining heat at the end of the routine (J.m-2) 93 96 REAL(wp), POINTER, DIMENSION(:) :: zf_tt ! Heat budget to determine melting or freezing(W.m-2) 94 INTEGER , POINTER, DIMENSION(:) :: icount ! number of layers vanished by melting95 97 96 98 REAL(wp), POINTER, DIMENSION(:) :: zdh_s_mel ! snow melt … … 100 102 REAL(wp), POINTER, DIMENSION(:,:) :: zdeltah 101 103 REAL(wp), POINTER, DIMENSION(:,:) :: zh_i ! ice layer thickness 104 INTEGER , POINTER, DIMENSION(:,:) :: icount ! number of layers vanished by melting 102 105 103 106 REAL(wp), POINTER, DIMENSION(:) :: zqh_i ! total ice heat content (J.m-2) 104 107 REAL(wp), POINTER, DIMENSION(:) :: zqh_s ! total snow heat content (J.m-2) 105 108 REAL(wp), POINTER, DIMENSION(:) :: zq_s ! total snow enthalpy (J.m-3) 109 REAL(wp), POINTER, DIMENSION(:) :: zsnw ! distribution of snow after wind blowing 106 110 107 111 REAL(wp) :: zswitch_sal … … 118 122 END SELECT 119 123 120 CALL wrk_alloc( jpij, z h_s, zqprec, zq_su, zq_bo, zf_tt, zq_rema)124 CALL wrk_alloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw ) 121 125 CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 122 CALL wrk_alloc( jpij, nlay_i +1, zdeltah, zh_i )123 CALL wrk_alloc( jpij, icount )124 126 CALL wrk_alloc( jpij, nlay_i, zdeltah, zh_i ) 127 CALL wrk_alloc( jpij, nlay_i, icount ) 128 125 129 dh_i_surf (:) = 0._wp ; dh_i_bott (:) = 0._wp ; dh_snowice(:) = 0._wp 126 130 dsm_i_se_1d(:) = 0._wp ; dsm_i_si_1d(:) = 0._wp 127 128 zqprec (:) = 0._wp ; zq_su (:) = 0._wp ; zq_bo (:) = 0._wp ; zf_tt(:) = 0._wp129 zq_rema (:) = 0._wp130 131 z h_s (:) = 0._wp132 zdh_s_pre(:) = 0._wp 133 zd h_s_mel(:) = 0._wp134 zdh_s_sub(:) = 0._wp135 zqh_s (:) = 0._wp 136 zqh_i (:) = 0._wp 137 138 zh_i (:,:) = 0._wp139 zdeltah (:,:) = 0._wp140 icount (:) = 0131 132 zqprec (:) = 0._wp ; zq_su (:) = 0._wp ; zq_bo (:) = 0._wp ; zf_tt(:) = 0._wp 133 zq_rema (:) = 0._wp ; zsnw (:) = 0._wp 134 zdh_s_mel(:) = 0._wp ; zdh_s_pre(:) = 0._wp ; zdh_s_sub(:) = 0._wp ; zqh_i(:) = 0._wp 135 zqh_s (:) = 0._wp ; zq_s (:) = 0._wp 136 137 zdeltah(:,:) = 0._wp ; zh_i(:,:) = 0._wp 138 icount (:,:) = 0 139 140 141 ! Initialize enthalpy at nlay_i+1 142 DO ji = kideb, kiut 143 q_i_1d(ji,nlay_i+1) = 0._wp 144 END DO 141 145 142 146 ! initialize layer thicknesses and enthalpies … … 155 159 ! 156 160 DO ji = kideb, kiut 157 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) )158 ztmelts = rswitch * rt0 + ( 1._wp - rswitch ) * rt0159 160 161 zfdum = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 161 162 zf_tt(ji) = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) 162 163 163 zq_su (ji) = MAX( 0._wp, zfdum * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - ztmelts) )164 zq_su (ji) = MAX( 0._wp, zfdum * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) ) 164 165 zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice ) 165 166 END DO … … 187 188 !------------------------------------------------------------! 188 189 ! 189 DO ji = kideb, kiut190 zh_s(ji) = ht_s_1d(ji) * r1_nlay_s191 END DO192 !193 190 DO jk = 1, nlay_s 194 191 DO ji = kideb, kiut 195 zqh_s(ji) = zqh_s(ji) + q_s_1d(ji,jk) * zh_s(ji)192 zqh_s(ji) = zqh_s(ji) + q_s_1d(ji,jk) * ht_s_1d(ji) * r1_nlay_s 196 193 END DO 197 194 END DO … … 222 219 ! Martin Vancoppenolle, December 2006 223 220 221 CALL lim_thd_snwblow( 1. - at_i_1d(kideb:kiut), zsnw(kideb:kiut) ) ! snow distribution over ice after wind blowing 222 223 zdeltah(:,:) = 0._wp 224 224 DO ji = kideb, kiut 225 225 !----------- … … 227 227 !----------- 228 228 ! thickness change 229 zcoeff = ( 1._wp - ( 1._wp - at_i_1d(ji) )**rn_betas ) / at_i_1d(ji) 230 zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice * r1_rhosn 231 ! enthalpy of the precip (>0, J.m-3) (tatm_ice is now in K) 232 zqprec (ji) = rhosn * ( cpic * ( rt0 - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus ) 229 zdh_s_pre(ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhosn / at_i_1d(ji) 230 ! enthalpy of the precip (>0, J.m-3) 231 zqprec (ji) = - qprec_ice_1d(ji) 233 232 IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp 234 233 ! heat flux from snow precip (>0, W.m-2) … … 236 235 ! mass flux, <0 237 236 wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice 238 ! update thickness239 ht_s_1d (ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_pre(ji) )240 237 241 238 !--------------------- … … 243 240 !--------------------- 244 241 ! thickness change 245 IF( zdh_s_pre(ji) > 0._wp ) THEN246 242 rswitch = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) ) 247 zd h_s_mel (ji) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 )248 zd h_s_mel (ji) = MAX( - zdh_s_pre(ji), zdh_s_mel(ji) ) ! bound melting243 zdeltah (ji,1) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 ) 244 zdeltah (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) ) ! bound melting 249 245 ! heat used to melt snow (W.m-2, >0) 250 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zd h_s_mel(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice246 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_rdtice 251 247 ! snow melting only = water into the ocean (then without snow precip), >0 252 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_mel(ji) * r1_rdtice 253 254 ! updates available heat + thickness 255 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdh_s_mel(ji) * zqprec(ji) ) 256 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_mel(ji) ) 257 zh_s (ji) = ht_s_1d(ji) * r1_nlay_s 258 259 ENDIF 260 END DO 261 262 ! If heat still available, then melt more snow 263 zdeltah(:,:) = 0._wp ! important 248 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice 249 ! updates available heat + precipitations after melting 250 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,1) * zqprec(ji) ) 251 zdh_s_pre (ji) = zdh_s_pre(ji) + zdeltah(ji,1) 252 253 ! update thickness 254 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_pre(ji) ) 255 END DO 256 257 ! If heat still available (zq_su > 0), then melt more snow 258 zdeltah(:,:) = 0._wp 264 259 DO jk = 1, nlay_s 265 260 DO ji = kideb, kiut … … 268 263 rswitch = rswitch * ( MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,jk) - epsi20 ) ) ) 269 264 zdeltah (ji,jk) = - rswitch * zq_su(ji) / MAX( q_s_1d(ji,jk), epsi20 ) 270 zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) ) ! bound melting265 zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - ht_s_1d(ji) ) ! bound melting 271 266 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) 272 267 ! heat used to melt snow(W.m-2, >0) … … 274 269 ! snow melting only = water into the ocean (then without snow precip) 275 270 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 276 277 271 ! updates available heat + thickness 278 272 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) ) 279 273 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdeltah(ji,jk) ) 280 281 274 END DO 282 275 END DO … … 286 279 !---------------------- 287 280 ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates 288 ! clem comment: not counted in mass exchange in limsbc since this is an exchange with atm. (not ocean)281 ! clem comment: not counted in mass/heat exchange in limsbc since this is an exchange with atm. (not ocean) 289 282 ! clem comment: ice should also sublimate 290 IF( lk_cpl ) THEN 291 ! coupled mode: sublimation already included in emp_ice (to do in limsbc_ice) 292 zdh_s_sub(:) = 0._wp 293 ELSE 294 ! forced mode: snow thickness change due to sublimation 295 DO ji = kideb, kiut 296 zdh_s_sub(ji) = MAX( - ht_s_1d(ji) , - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice ) 297 ! Heat flux by sublimation [W.m-2], < 0 298 ! sublimate first snow that had fallen, then pre-existing snow 299 zcoeff = ( MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) ) * zqprec(ji) + & 300 & ( zdh_s_sub(ji) - MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) ) ) * q_s_1d(ji,1) ) & 301 & * a_i_1d(ji) * r1_rdtice 302 hfx_sub_1d(ji) = hfx_sub_1d(ji) + zcoeff 303 ! Mass flux by sublimation 304 wfx_sub_1d(ji) = wfx_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice 305 ! new snow thickness 306 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) ) 307 END DO 308 ENDIF 309 283 zdeltah(:,:) = 0._wp 284 ! coupled mode: sublimation is set to 0 (evap_ice = 0) until further notice 285 ! forced mode: snow thickness change due to sublimation 286 DO ji = kideb, kiut 287 zdh_s_sub(ji) = MAX( - ht_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice ) 288 ! Heat flux by sublimation [W.m-2], < 0 289 ! sublimate first snow that had fallen, then pre-existing snow 290 zdeltah(ji,1) = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) ) 291 hfx_sub_1d(ji) = hfx_sub_1d(ji) + ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * q_s_1d(ji,1) & 292 & ) * a_i_1d(ji) * r1_rdtice 293 ! Mass flux by sublimation 294 wfx_sub_1d(ji) = wfx_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice 295 ! new snow thickness 296 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) ) 297 ! update precipitations after sublimation and correct sublimation 298 zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1) 299 zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1) 300 END DO 301 310 302 ! --- Update snow diags --- ! 311 303 DO ji = kideb, kiut 312 dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji) 313 zh_s(ji) = ht_s_1d(ji) * r1_nlay_s 314 END DO ! ji 304 dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji) 305 END DO 315 306 316 307 !------------------------------------------- … … 323 314 rswitch = MAX( 0._wp , SIGN( 1._wp, ht_s_1d(ji) - epsi20 ) ) 324 315 q_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) * & 325 & ( ( MAX( 0._wp, dh_s_tot(ji) )) * zqprec(ji) + &326 & ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_1d(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) )316 & ( ( zdh_s_pre(ji) ) * zqprec(ji) + & 317 & ( ht_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) ) 327 318 zq_s(ji) = zq_s(ji) + q_s_1d(ji,jk) 328 319 END DO … … 334 325 zdeltah(:,:) = 0._wp ! important 335 326 DO jk = 1, nlay_i 336 DO ji = kideb, kiut 337 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of layer k [J/kg, <0] 338 339 ztmelts = - tmut * s_i_1d(ji,jk) + rt0 ! Melting point of layer k [K] 340 341 zEw = rcp * ( ztmelts - rt0 ) ! Specific enthalpy of resulting meltwater [J/kg, <0] 342 343 zdE = zEi - zEw ! Specific enthalpy difference < 0 344 345 zfmdt = - zq_su(ji) / zdE ! Mass flux to the ocean [kg/m2, >0] 346 347 zdeltah(ji,jk) = - zfmdt * r1_rhoic ! Melt of layer jk [m, <0] 348 349 zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk) , - zh_i(ji,jk) ) ) ! Melt of layer jk cannot exceed the layer thickness [m, <0] 350 351 zq_su(ji) = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat 352 353 dh_i_surf(ji) = dh_i_surf(ji) + zdeltah(ji,jk) ! Cumulate surface melt 354 355 zfmdt = - rhoic * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0] 356 357 zQm = zfmdt * zEw ! Energy of the melt water sent to the ocean [J/m2, <0] 358 359 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 360 sfx_sum_1d(ji) = sfx_sum_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice 361 362 ! Contribution to heat flux [W.m-2], < 0 363 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 364 365 ! Total heat flux used in this process [W.m-2], > 0 366 hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 367 368 ! Contribution to mass flux 369 wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 370 327 DO ji = kideb, kiut 328 ztmelts = - tmut * s_i_1d(ji,jk) + rt0 ! Melting point of layer k [K] 329 330 IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting 331 332 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of layer k [J/kg, <0] 333 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) 334 ! set up at 0 since no energy is needed to melt water...(it is already melted) 335 zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing 336 ! this should normally not happen, but sometimes, heat diffusion leads to this 337 zfmdt = - zdeltah(ji,jk) * rhoic ! Mass flux x time step > 0 338 339 dh_i_surf(ji) = dh_i_surf(ji) + zdeltah(ji,jk) ! Cumulate surface melt 340 341 zfmdt = - rhoic * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0] 342 343 ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean) 344 hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice 345 346 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 347 sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 348 349 ! Contribution to mass flux 350 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 351 352 ELSE !!! Surface melting 353 354 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of layer k [J/kg, <0] 355 zEw = rcp * ( ztmelts - rt0 ) ! Specific enthalpy of resulting meltwater [J/kg, <0] 356 zdE = zEi - zEw ! Specific enthalpy difference < 0 357 358 zfmdt = - zq_su(ji) / zdE ! Mass flux to the ocean [kg/m2, >0] 359 360 zdeltah(ji,jk) = - zfmdt * r1_rhoic ! Melt of layer jk [m, <0] 361 362 zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk) , - zh_i(ji,jk) ) ) ! Melt of layer jk cannot exceed the layer thickness [m, <0] 363 364 zq_su(ji) = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat 365 366 dh_i_surf(ji) = dh_i_surf(ji) + zdeltah(ji,jk) ! Cumulate surface melt 367 368 zfmdt = - rhoic * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0] 369 370 zQm = zfmdt * zEw ! Energy of the melt water sent to the ocean [J/m2, <0] 371 372 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 373 sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 374 375 ! Contribution to heat flux [W.m-2], < 0 376 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 377 378 ! Total heat flux used in this process [W.m-2], > 0 379 hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 380 381 ! Contribution to mass flux 382 wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 383 384 END IF 371 385 ! record which layers have disappeared (for bottom melting) 372 386 ! => icount=0 : no layer has vanished 373 387 ! => icount=5 : 5 layers have vanished 374 rswitch = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) )375 icount(ji ) = icount(ji) +NINT( rswitch )376 zh_i(ji,jk) = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) )388 rswitch = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) ) 389 icount(ji,jk) = NINT( rswitch ) 390 zh_i(ji,jk) = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) ) 377 391 378 392 ! update heat content (J.m-2) and layer thickness … … 405 419 ! -> need for an iterative procedure, which converges quickly 406 420 407 IF ( nn_icesal == 2 ) THEN 408 num_iter_max = 5 409 ELSE 410 num_iter_max = 1 411 ENDIF 412 413 ! Just to be sure that enthalpy at nlay_i+1 is null 414 DO ji = kideb, kiut 415 q_i_1d(ji,nlay_i+1) = 0._wp 416 END DO 421 num_iter_max = 1 422 IF( nn_icesal == 2 ) num_iter_max = 5 417 423 418 424 ! Iterative procedure … … 483 489 484 490 ! Contribution to salt flux, <0 485 sfx_bog_1d(ji) = sfx_bog_1d(ji) + s_i_new(ji) * a_i_1d(ji) * zfmdt* r1_rdtice491 sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * s_i_new(ji) * r1_rdtice 486 492 487 493 ! Contribution to mass flux, <0 … … 500 506 DO jk = nlay_i, 1, -1 501 507 DO ji = kideb, kiut 502 IF( zf_tt(ji) > = 0._wp .AND. jk > icount(ji) ) THEN ! do not calculate where layer has already disappeared by surface melting508 IF( zf_tt(ji) > 0._wp .AND. jk > icount(ji,jk) ) THEN ! do not calculate where layer has already disappeared by surface melting 503 509 504 510 ztmelts = - tmut * s_i_1d(ji,jk) + rt0 ! Melting point of layer jk (K) … … 507 513 508 514 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) 509 510 !!zEw = rcp * ( t_i_1d(ji,jk) - rt0 ) ! Specific enthalpy of meltwater at T = t_i_1d (J/kg, <0)511 512 515 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) 513 516 ! set up at 0 since no energy is needed to melt water...(it is already melted) 514 515 zdeltah (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing 516 ! this should normally not happen, but sometimes, heat diffusion leads to this 517 zdeltah (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing 518 ! this should normally not happen, but sometimes, heat diffusion leads to this 517 519 518 520 dh_i_bott (ji) = dh_i_bott(ji) + zdeltah(ji,jk) 519 521 520 zfmdt = - zdeltah(ji,jk) * rhoic 522 zfmdt = - zdeltah(ji,jk) * rhoic ! Mass flux x time step > 0 521 523 522 524 ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean) … … 524 526 525 527 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 526 sfx_res_1d(ji) = sfx_res_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic* r1_rdtice528 sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 527 529 528 530 ! Contribution to mass flux … … 535 537 ELSE !!! Basal melting 536 538 537 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) 538 539 zEw = rcp * ( ztmelts - rt0 ) ! Specific enthalpy of meltwater (J/kg, <0) 540 541 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) 542 543 zfmdt = - zq_bo(ji) / zdE ! Mass flux x time step (kg/m2, >0) 544 545 zdeltah(ji,jk) = - zfmdt * r1_rhoic ! Gross thickness change 546 547 zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change 539 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) 540 zEw = rcp * ( ztmelts - rt0 ) ! Specific enthalpy of meltwater (J/kg, <0) 541 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) 542 543 zfmdt = - zq_bo(ji) / zdE ! Mass flux x time step (kg/m2, >0) 544 545 zdeltah(ji,jk) = - zfmdt * r1_rhoic ! Gross thickness change 546 547 zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change 548 548 549 zq_bo(ji) 550 551 dh_i_bott(ji) 552 553 zfmdt 554 555 zQm 549 zq_bo(ji) = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat. MAX is necessary for roundup errors 550 551 dh_i_bott(ji) = dh_i_bott(ji) + zdeltah(ji,jk) ! Update basal melt 552 553 zfmdt = - zdeltah(ji,jk) * rhoic ! Mass flux x time step > 0 554 555 zQm = zfmdt * zEw ! Heat exchanged with ocean 556 556 557 557 ! Contribution to heat flux to the ocean [W.m-2], <0 558 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice558 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 559 559 560 560 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 561 sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic* r1_rdtice561 sfx_bom_1d(ji) = sfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 562 562 563 563 ! Total heat flux used in this process [W.m-2], >0 564 hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice564 hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 565 565 566 566 ! Contribution to mass flux 567 wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice567 wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 568 568 569 569 ! update heat content (J.m-2) and layer thickness … … 595 595 zdeltah (ji,1) = - rswitch * zq_rema(ji) / MAX( q_s_1d(ji,1), epsi20 ) 596 596 zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_1d(ji) ) ) ! bound melting 597 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,1)598 597 dh_s_tot (ji) = dh_s_tot(ji) + zdeltah(ji,1) 599 598 ht_s_1d (ji) = ht_s_1d(ji) + zdeltah(ji,1) … … 622 621 dh_snowice(ji) = MAX( 0._wp , ( rhosn * ht_s_1d(ji) + (rhoic-rau0) * ht_i_1d(ji) ) / ( rhosn+rau0-rhoic ) ) 623 622 624 ht_i_1d(ji) 625 ht_s_1d(ji) 623 ht_i_1d(ji) = ht_i_1d(ji) + dh_snowice(ji) 624 ht_s_1d(ji) = ht_s_1d(ji) - dh_snowice(ji) 626 625 627 626 ! Salinity of snow ice … … 669 668 ! Update temperature, energy 670 669 !------------------------------------------- 671 !clem bug: we should take snow into account here672 670 DO ji = kideb, kiut 673 671 rswitch = 1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) … … 688 686 WHERE( ht_i_1d == 0._wp ) a_i_1d = 0._wp 689 687 690 CALL wrk_dealloc( jpij, z h_s, zqprec, zq_su, zq_bo, zf_tt, zq_rema)688 CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw ) 691 689 CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 692 CALL wrk_dealloc( jpij, nlay_i +1, zdeltah, zh_i )693 CALL wrk_dealloc( jpij, icount )690 CALL wrk_dealloc( jpij, nlay_i, zdeltah, zh_i ) 691 CALL wrk_dealloc( jpij, nlay_i, icount ) 694 692 ! 695 693 ! 696 694 END SUBROUTINE lim_thd_dh 695 696 697 !!-------------------------------------------------------------------------- 698 !! INTERFACE lim_thd_snwblow 699 !! ** Purpose : Compute distribution of precip over the ice 700 !!-------------------------------------------------------------------------- 701 SUBROUTINE lim_thd_snwblow_2d( pin, pout ) 702 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pin ! previous fraction lead ( pfrld or (1. - a_i_b) ) 703 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pout 704 pout = ( 1._wp - ( pin )**rn_betas ) 705 END SUBROUTINE lim_thd_snwblow_2d 706 707 SUBROUTINE lim_thd_snwblow_1d( pin, pout ) 708 REAL(wp), DIMENSION(:), INTENT(in ) :: pin 709 REAL(wp), DIMENSION(:), INTENT(inout) :: pout 710 pout = ( 1._wp - ( pin )**rn_betas ) 711 END SUBROUTINE lim_thd_snwblow_1d 712 697 713 698 714 #else
Note: See TracChangeset
for help on using the changeset viewer.
