- Timestamp:
- 2015-06-22T16:40:58+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/restart_datestamp/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r5420 r5462 20 20 USE sbc_oce ! Surface boundary condition: ocean fields 21 21 USE ice ! LIM variables 22 USE par_ice ! LIM parameters23 22 USE thd_ice ! LIM thermodynamics 24 23 USE in_out_manager ! I/O manager … … 30 29 PRIVATE 31 30 32 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/sbccpl and here 33 34 INTERFACE lim_thd_snwblow 35 MODULE PROCEDURE lim_thd_snwblow_1d, lim_thd_snwblow_2d 36 END INTERFACE 33 37 34 38 !!---------------------------------------------------------------------- … … 70 74 71 75 REAL(wp) :: ztmelts ! local scalar 72 REAL(wp) :: z dh, zfdum !76 REAL(wp) :: zfdum 73 77 REAL(wp) :: zfracs ! fractionation coefficient for bottom salt entrapment 74 REAL(wp) :: zcoeff ! dummy argument for snowfall partitioning over ice and leads 75 REAL(wp) :: zs_snic ! snow-ice salinity 78 REAL(wp) :: zs_snic ! snow-ice salinity 76 79 REAL(wp) :: zswi1 ! switch for computation of bottom salinity 77 80 REAL(wp) :: zswi12 ! switch for computation of bottom salinity … … 87 90 REAL(wp) :: zsstK ! SST in Kelvin 88 91 89 REAL(wp), POINTER, DIMENSION(:) :: zh_s ! snow layer thickness90 92 REAL(wp), POINTER, DIMENSION(:) :: zqprec ! energy of fallen snow (J.m-3) 91 93 REAL(wp), POINTER, DIMENSION(:) :: zq_su ! heat for surface ablation (J.m-2) 92 94 REAL(wp), POINTER, DIMENSION(:) :: zq_bo ! heat for bottom ablation (J.m-2) 93 REAL(wp), POINTER, DIMENSION(:) :: zq_1cat ! corrected heat in case 1-cat and hmelt>15cm (J.m-2)94 95 REAL(wp), POINTER, DIMENSION(:) :: zq_rema ! remaining heat at the end of the routine (J.m-2) 95 REAL(wp), POINTER, DIMENSION(:) :: zf_tt ! Heat budget to determine melting or freezing(W.m-2) 96 INTEGER , POINTER, DIMENSION(:) :: icount ! number of layers vanished by melting 96 REAL(wp), POINTER, DIMENSION(:) :: zf_tt ! Heat budget to determine melting or freezing(W.m-2) 97 97 98 98 REAL(wp), POINTER, DIMENSION(:) :: zdh_s_mel ! snow melt … … 102 102 REAL(wp), POINTER, DIMENSION(:,:) :: zdeltah 103 103 REAL(wp), POINTER, DIMENSION(:,:) :: zh_i ! ice layer thickness 104 INTEGER , POINTER, DIMENSION(:,:) :: icount ! number of layers vanished by melting 104 105 105 106 REAL(wp), POINTER, DIMENSION(:) :: zqh_i ! total ice heat content (J.m-2) 106 107 REAL(wp), POINTER, DIMENSION(:) :: zqh_s ! total snow heat content (J.m-2) 107 108 REAL(wp), POINTER, DIMENSION(:) :: zq_s ! total snow enthalpy (J.m-3) 108 109 ! mass and salt flux (clem) 110 REAL(wp) :: z dvres, zswitch_sal109 REAL(wp), POINTER, DIMENSION(:) :: zsnw ! distribution of snow after wind blowing 110 111 REAL(wp) :: zswitch_sal 111 112 112 113 ! Heat conservation … … 115 116 !!------------------------------------------------------------------ 116 117 117 ! Discriminate between varying salinity (n um_sal=2) and prescribed cases (other values)118 SELECT CASE( n um_sal ) ! varying salinity or not118 ! Discriminate between varying salinity (nn_icesal=2) and prescribed cases (other values) 119 SELECT CASE( nn_icesal ) ! varying salinity or not 119 120 CASE( 1, 3, 4 ) ; zswitch_sal = 0 ! prescribed salinity profile 120 121 CASE( 2 ) ; zswitch_sal = 1 ! varying salinity profile 121 122 END SELECT 122 123 123 CALL wrk_alloc( jpij, z h_s, zqprec, zq_su, zq_bo, zf_tt, zq_1cat, zq_rema )124 CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s )125 CALL wrk_alloc( jpij, nlay_i +1, zdeltah, zh_i )126 CALL wrk_alloc( jpij, icount )124 CALL wrk_alloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema ) 125 CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s, zsnw ) 126 CALL wrk_alloc( jpij, nlay_i, zdeltah, zh_i ) 127 CALL wrk_alloc( jpij, nlay_i, icount ) 127 128 128 129 dh_i_surf (:) = 0._wp ; dh_i_bott (:) = 0._wp ; dh_snowice(:) = 0._wp … … 130 131 131 132 zqprec (:) = 0._wp ; zq_su (:) = 0._wp ; zq_bo (:) = 0._wp ; zf_tt (:) = 0._wp 132 zq_1cat(:) = 0._wp ; zq_rema(:) = 0._wp 133 134 zh_s (:) = 0._wp 133 zq_rema(:) = 0._wp 134 135 135 zdh_s_pre(:) = 0._wp 136 136 zdh_s_mel(:) = 0._wp … … 141 141 zh_i (:,:) = 0._wp 142 142 zdeltah (:,:) = 0._wp 143 icount (:) = 0 143 icount (:,:) = 0 144 145 ! Initialize enthalpy at nlay_i+1 146 DO ji = kideb, kiut 147 q_i_1d(ji,nlay_i+1) = 0._wp 148 END DO 144 149 145 150 ! initialize layer thicknesses and enthalpies … … 148 153 DO jk = 1, nlay_i 149 154 DO ji = kideb, kiut 150 h_i_old (ji,jk) = ht_i_1d(ji) / REAL( nlay_i )155 h_i_old (ji,jk) = ht_i_1d(ji) * r1_nlay_i 151 156 qh_i_old(ji,jk) = q_i_1d(ji,jk) * h_i_old(ji,jk) 152 157 ENDDO … … 158 163 ! 159 164 DO ji = kideb, kiut 160 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_1d(ji) ) )161 ztmelts = rswitch * rtt + ( 1._wp - rswitch ) * rtt162 163 165 zfdum = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 164 166 zf_tt(ji) = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) 165 167 166 zq_su (ji) = MAX( 0._wp, zfdum * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - ztmelts) )168 zq_su (ji) = MAX( 0._wp, zfdum * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) ) 167 169 zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice ) 168 170 END DO … … 174 176 !------------------------------------------------------------------------------! 175 177 DO ji = kideb, kiut 176 IF( t_s_1d(ji,1) > rt t) THEN !!! Internal melting178 IF( t_s_1d(ji,1) > rt0 ) THEN !!! Internal melting 177 179 ! Contribution to heat flux to the ocean [W.m-2], < 0 178 180 hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_1d(ji,1) * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice … … 182 184 ht_s_1d(ji) = 0._wp 183 185 q_s_1d (ji,1) = 0._wp 184 t_s_1d (ji,1) = rt t186 t_s_1d (ji,1) = rt0 185 187 END IF 186 188 END DO … … 190 192 !------------------------------------------------------------! 191 193 ! 192 DO ji = kideb, kiut193 zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s )194 END DO195 !196 194 DO jk = 1, nlay_s 197 195 DO ji = kideb, kiut 198 zqh_s(ji) = zqh_s(ji) + q_s_1d(ji,jk) * zh_s(ji)196 zqh_s(ji) = zqh_s(ji) + q_s_1d(ji,jk) * ht_s_1d(ji) * r1_nlay_s 199 197 END DO 200 198 END DO … … 202 200 DO jk = 1, nlay_i 203 201 DO ji = kideb, kiut 204 zh_i(ji,jk) = ht_i_1d(ji) / REAL( nlay_i )202 zh_i(ji,jk) = ht_i_1d(ji) * r1_nlay_i 205 203 zqh_i(ji) = zqh_i(ji) + q_i_1d(ji,jk) * zh_i(ji,jk) 206 204 END DO … … 225 223 ! Martin Vancoppenolle, December 2006 226 224 225 zdeltah(:,:) = 0._wp 226 CALL lim_thd_snwblow( 1. - at_i_1d, zsnw ) ! snow distribution over ice after wind blowing 227 227 DO ji = kideb, kiut 228 228 !----------- … … 230 230 !----------- 231 231 ! thickness change 232 zcoeff = ( 1._wp - ( 1._wp - at_i_1d(ji) )**betas ) / at_i_1d(ji) 233 zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice / rhosn 234 ! enthalpy of the precip (>0, J.m-3) (tatm_ice is now in K) 235 zqprec (ji) = rhosn * ( cpic * ( rtt - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus ) 232 zdh_s_pre(ji) = zsnw(ji) * sprecip_1d(ji) * rdt_ice * r1_rhosn / at_i_1d(ji) 233 ! enthalpy of the precip (>0, J.m-3) 234 zqprec (ji) = - qprec_ice_1d(ji) 236 235 IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp 237 236 ! heat flux from snow precip (>0, W.m-2) … … 239 238 ! mass flux, <0 240 239 wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_pre(ji) * r1_rdtice 241 ! update thickness242 ht_s_1d (ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_pre(ji) )243 240 244 241 !--------------------- … … 246 243 !--------------------- 247 244 ! thickness change 248 IF( zdh_s_pre(ji) > 0._wp ) THEN 249 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zqprec(ji) + epsi20 ) ) 250 zdh_s_mel (ji) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 ) 251 zdh_s_mel (ji) = MAX( - zdh_s_pre(ji), zdh_s_mel(ji) ) ! bound melting 245 rswitch = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) ) 246 zdeltah (ji,1) = - rswitch * zq_su(ji) / MAX( zqprec(ji) , epsi20 ) 247 zdeltah (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) ) ! bound melting 252 248 ! heat used to melt snow (W.m-2, >0) 253 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zd h_s_mel(ji) * a_i_1d(ji) * zqprec(ji) * r1_rdtice249 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_rdtice 254 250 ! snow melting only = water into the ocean (then without snow precip), >0 255 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_mel(ji) * r1_rdtice 256 257 ! updates available heat + thickness 258 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdh_s_mel(ji) * zqprec(ji) ) 259 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_mel(ji) ) 260 zh_s (ji) = ht_s_1d(ji) / REAL( nlay_s ) 261 262 ENDIF 263 END DO 264 265 ! If heat still available, then melt more snow 266 zdeltah(:,:) = 0._wp ! important 251 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice 252 ! updates available heat + precipitations after melting 253 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,1) * zqprec(ji) ) 254 zdh_s_pre (ji) = zdh_s_pre(ji) + zdeltah(ji,1) 255 256 ! update thickness 257 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_pre(ji) ) 258 END DO 259 260 ! If heat still available (zq_su > 0), then melt more snow 261 zdeltah(:,:) = 0._wp 267 262 DO jk = 1, nlay_s 268 263 DO ji = kideb, kiut 269 264 ! thickness change 270 265 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) ) 271 rswitch = rswitch * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, - q_s_1d(ji,jk) +epsi20 ) ) )266 rswitch = rswitch * ( MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,jk) - epsi20 ) ) ) 272 267 zdeltah (ji,jk) = - rswitch * zq_su(ji) / MAX( q_s_1d(ji,jk), epsi20 ) 273 zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) ) ! bound melting268 zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - ht_s_1d(ji) ) ! bound melting 274 269 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) 275 270 ! heat used to melt snow(W.m-2, >0) … … 277 272 ! snow melting only = water into the ocean (then without snow precip) 278 273 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 279 280 274 ! updates available heat + thickness 281 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) )275 zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) ) 282 276 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdeltah(ji,jk) ) 283 284 277 END DO 285 278 END DO … … 289 282 !---------------------- 290 283 ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates 291 ! clem comment: not counted in mass exchange in limsbc since this is an exchange with atm. (not ocean)284 ! clem comment: not counted in mass/heat exchange in limsbc since this is an exchange with atm. (not ocean) 292 285 ! clem comment: ice should also sublimate 293 IF( lk_cpl ) THEN 294 ! coupled mode: sublimation already included in emp_ice (to do in limsbc_ice) 295 zdh_s_sub(:) = 0._wp 296 ELSE 297 ! forced mode: snow thickness change due to sublimation 298 DO ji = kideb, kiut 299 zdh_s_sub(ji) = MAX( - ht_s_1d(ji) , - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice ) 300 ! Heat flux by sublimation [W.m-2], < 0 301 ! sublimate first snow that had fallen, then pre-existing snow 302 zcoeff = ( MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) ) * zqprec(ji) + & 303 & ( 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) ) & 304 & * a_i_1d(ji) * r1_rdtice 305 hfx_sub_1d(ji) = hfx_sub_1d(ji) + zcoeff 306 ! Mass flux by sublimation 307 wfx_sub_1d(ji) = wfx_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice 308 ! new snow thickness 309 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) ) 310 END DO 311 ENDIF 312 286 zdeltah(:,:) = 0._wp 287 ! coupled mode: sublimation is set to 0 (evap_ice = 0) until further notice 288 ! forced mode: snow thickness change due to sublimation 289 DO ji = kideb, kiut 290 zdh_s_sub(ji) = MAX( - ht_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice ) 291 ! Heat flux by sublimation [W.m-2], < 0 292 ! sublimate first snow that had fallen, then pre-existing snow 293 zdeltah(ji,1) = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) ) 294 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) & 295 & ) * a_i_1d(ji) * r1_rdtice 296 ! Mass flux by sublimation 297 wfx_sub_1d(ji) = wfx_sub_1d(ji) - rhosn * a_i_1d(ji) * zdh_s_sub(ji) * r1_rdtice 298 ! new snow thickness 299 ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_sub(ji) ) 300 ! update precipitations after sublimation and correct sublimation 301 zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1) 302 zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1) 303 END DO 304 313 305 ! --- Update snow diags --- ! 314 306 DO ji = kideb, kiut 315 dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji) 316 zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s ) 317 END DO ! ji 307 dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji) 308 END DO 318 309 319 310 !------------------------------------------- … … 324 315 DO jk = 1, nlay_s 325 316 DO ji = kideb,kiut 326 rswitch = MAX( 0._wp , SIGN( 1._wp, - ht_s_1d(ji) +epsi20 ) )327 q_s_1d(ji,jk) = ( 1._wp - rswitch ) / MAX( ht_s_1d(ji), epsi20 ) *&328 & ( ( MAX( 0._wp, dh_s_tot(ji) )) * zqprec(ji) + &329 & ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_1d(ji) ) * rhosn * ( cpic * ( rtt- t_s_1d(ji,jk) ) + lfus ) )317 rswitch = MAX( 0._wp , SIGN( 1._wp, ht_s_1d(ji) - epsi20 ) ) 318 q_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) * & 319 & ( ( zdh_s_pre(ji) ) * zqprec(ji) + & 320 & ( ht_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) ) 330 321 zq_s(ji) = zq_s(ji) + q_s_1d(ji,jk) 331 322 END DO … … 337 328 zdeltah(:,:) = 0._wp ! important 338 329 DO jk = 1, nlay_i 339 DO ji = kideb, kiut 340 zEi = - q_i_1d(ji,jk) / rhoic ! Specific enthalpy of layer k [J/kg, <0] 341 342 ztmelts = - tmut * s_i_1d(ji,jk) + rtt ! Melting point of layer k [K] 343 344 zEw = rcp * ( ztmelts - rt0 ) ! Specific enthalpy of resulting meltwater [J/kg, <0] 345 346 zdE = zEi - zEw ! Specific enthalpy difference < 0 347 348 zfmdt = - zq_su(ji) / zdE ! Mass flux to the ocean [kg/m2, >0] 349 350 zdeltah(ji,jk) = - zfmdt / rhoic ! Melt of layer jk [m, <0] 351 352 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] 353 354 zq_su(ji) = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat 355 356 dh_i_surf(ji) = dh_i_surf(ji) + zdeltah(ji,jk) ! Cumulate surface melt 357 358 zfmdt = - rhoic * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0] 359 360 zQm = zfmdt * zEw ! Energy of the melt water sent to the ocean [J/m2, <0] 361 362 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 363 sfx_sum_1d(ji) = sfx_sum_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice 364 365 ! Contribution to heat flux [W.m-2], < 0 366 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 367 368 ! Total heat flux used in this process [W.m-2], > 0 369 hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 370 371 ! Contribution to mass flux 372 wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 373 330 DO ji = kideb, kiut 331 ztmelts = - tmut * s_i_1d(ji,jk) + rt0 ! Melting point of layer k [K] 332 333 IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting 334 335 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of layer k [J/kg, <0] 336 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) 337 ! set up at 0 since no energy is needed to melt water...(it is already melted) 338 zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing 339 ! this should normally not happen, but sometimes, heat diffusion leads to this 340 zfmdt = - zdeltah(ji,jk) * rhoic ! Mass flux x time step > 0 341 342 dh_i_surf(ji) = dh_i_surf(ji) + zdeltah(ji,jk) ! Cumulate surface melt 343 344 zfmdt = - rhoic * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0] 345 346 ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean) 347 hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_rdtice 348 349 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 350 sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 351 352 ! Contribution to mass flux 353 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 354 355 ELSE !!! Surface melting 356 357 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of layer k [J/kg, <0] 358 zEw = rcp * ( ztmelts - rt0 ) ! Specific enthalpy of resulting meltwater [J/kg, <0] 359 zdE = zEi - zEw ! Specific enthalpy difference < 0 360 361 zfmdt = - zq_su(ji) / zdE ! Mass flux to the ocean [kg/m2, >0] 362 363 zdeltah(ji,jk) = - zfmdt * r1_rhoic ! Melt of layer jk [m, <0] 364 365 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] 366 367 zq_su(ji) = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat 368 369 dh_i_surf(ji) = dh_i_surf(ji) + zdeltah(ji,jk) ! Cumulate surface melt 370 371 zfmdt = - rhoic * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0] 372 373 zQm = zfmdt * zEw ! Energy of the melt water sent to the ocean [J/m2, <0] 374 375 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 376 sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 377 378 ! Contribution to heat flux [W.m-2], < 0 379 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 380 381 ! Total heat flux used in this process [W.m-2], > 0 382 hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 383 384 ! Contribution to mass flux 385 wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 386 387 END IF 374 388 ! record which layers have disappeared (for bottom melting) 375 389 ! => icount=0 : no layer has vanished 376 390 ! => icount=5 : 5 layers have vanished 377 rswitch = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) )378 icount(ji ) = icount(ji) +NINT( rswitch )379 zh_i(ji,jk) = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) )391 rswitch = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) ) 392 icount(ji,jk) = NINT( rswitch ) 393 zh_i(ji,jk) = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) ) 380 394 381 395 ! update heat content (J.m-2) and layer thickness … … 408 422 ! -> need for an iterative procedure, which converges quickly 409 423 410 IF ( num_sal == 2 ) THEN 411 num_iter_max = 5 412 ELSE 413 num_iter_max = 1 414 ENDIF 415 416 !clem debug. Just to be sure that enthalpy at nlay_i+1 is null 417 DO ji = kideb, kiut 418 q_i_1d(ji,nlay_i+1) = 0._wp 419 END DO 424 num_iter_max = 1 425 IF( nn_icesal == 2 ) num_iter_max = 5 420 426 421 427 ! Iterative procedure … … 440 446 + ( 1. - zswitch_sal ) * sm_i_1d(ji) 441 447 ! New ice growth 442 ztmelts = - tmut * s_i_new(ji) + rt t! New ice melting point (K)448 ztmelts = - tmut * s_i_new(ji) + rt0 ! New ice melting point (K) 443 449 444 450 zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 445 451 446 452 zEi = cpic * ( zt_i_new - ztmelts ) & ! Specific enthalpy of forming ice (J/kg, <0) 447 & - lfus * ( 1.0 - ( ztmelts - rt t ) / ( zt_i_new - rtt) ) &448 & + rcp * ( ztmelts-rt t)453 & - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) ) & 454 & + rcp * ( ztmelts-rt0 ) 449 455 450 456 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0) … … 456 462 q_i_1d(ji,nlay_i+1) = -zEi * rhoic ! New ice energy of melting (J/m3, >0) 457 463 458 ENDIF ! fc_bo_i459 END DO ! ji460 END DO ! iter464 ENDIF 465 END DO 466 END DO 461 467 462 468 ! Contribution to Energy and Salt Fluxes … … 467 473 zfmdt = - rhoic * dh_i_bott(ji) ! Mass flux x time step (kg/m2, < 0) 468 474 469 ztmelts = - tmut * s_i_new(ji) + rt t! New ice melting point (K)475 ztmelts = - tmut * s_i_new(ji) + rt0 ! New ice melting point (K) 470 476 471 477 zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 472 478 473 479 zEi = cpic * ( zt_i_new - ztmelts ) & ! Specific enthalpy of forming ice (J/kg, <0) 474 & - lfus * ( 1.0 - ( ztmelts - rt t ) / ( zt_i_new - rtt) ) &475 & + rcp * ( ztmelts-rt t)480 & - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) ) & 481 & + rcp * ( ztmelts-rt0 ) 476 482 477 483 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0) … … 486 492 487 493 ! Contribution to salt flux, <0 488 sfx_bog_1d(ji) = sfx_bog_1d(ji) + s_i_new(ji) * a_i_1d(ji) * zfmdt* r1_rdtice494 sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoic * a_i_1d(ji) * dh_i_bott(ji) * s_i_new(ji) * r1_rdtice 489 495 490 496 ! Contribution to mass flux, <0 … … 503 509 DO jk = nlay_i, 1, -1 504 510 DO ji = kideb, kiut 505 IF( zf_tt(ji) > = 0._wp .AND. jk > icount(ji) ) THEN ! do not calculate where layer has already disappeared fromsurface melting506 507 ztmelts = - tmut * s_i_1d(ji,jk) + rt t! Melting point of layer jk (K)511 IF( zf_tt(ji) > 0._wp .AND. jk > icount(ji,jk) ) THEN ! do not calculate where layer has already disappeared by surface melting 512 513 ztmelts = - tmut * s_i_1d(ji,jk) + rt0 ! Melting point of layer jk (K) 508 514 509 515 IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting 510 516 511 zEi = - q_i_1d(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0) 512 513 !!zEw = rcp * ( t_i_1d(ji,jk) - rtt ) ! Specific enthalpy of meltwater at T = t_i_1d (J/kg, <0) 514 517 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) 515 518 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) 516 519 ! set up at 0 since no energy is needed to melt water...(it is already melted) 517 518 zdeltah (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing 519 ! this should normally not happen, but sometimes, heat diffusion leads to this 520 zdeltah (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing 521 ! this should normally not happen, but sometimes, heat diffusion leads to this 520 522 521 523 dh_i_bott (ji) = dh_i_bott(ji) + zdeltah(ji,jk) 522 524 523 zfmdt = - zdeltah(ji,jk) * rhoic 525 zfmdt = - zdeltah(ji,jk) * rhoic ! Mass flux x time step > 0 524 526 525 527 ! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean) … … 527 529 528 530 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 529 sfx_res_1d(ji) = sfx_res_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic* r1_rdtice531 sfx_res_1d(ji) = sfx_res_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 530 532 531 533 ! Contribution to mass flux … … 538 540 ELSE !!! Basal melting 539 541 540 zEi = - q_i_1d(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0) 541 542 zEw = rcp * ( ztmelts - rtt )! Specific enthalpy of meltwater (J/kg, <0) 543 544 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) 545 546 zfmdt = - zq_bo(ji) / zdE ! Mass flux x time step (kg/m2, >0) 547 548 zdeltah(ji,jk) = - zfmdt / rhoic ! Gross thickness change 549 550 zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change 542 zEi = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0) 543 zEw = rcp * ( ztmelts - rt0 ) ! Specific enthalpy of meltwater (J/kg, <0) 544 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) 545 546 zfmdt = - zq_bo(ji) / zdE ! Mass flux x time step (kg/m2, >0) 547 548 zdeltah(ji,jk) = - zfmdt * r1_rhoic ! Gross thickness change 549 550 zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change 551 551 552 zq_bo(ji) 553 554 dh_i_bott(ji) 555 556 zfmdt 557 558 zQm 552 zq_bo(ji) = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat. MAX is necessary for roundup errors 553 554 dh_i_bott(ji) = dh_i_bott(ji) + zdeltah(ji,jk) ! Update basal melt 555 556 zfmdt = - zdeltah(ji,jk) * rhoic ! Mass flux x time step > 0 557 558 zQm = zfmdt * zEw ! Heat exchanged with ocean 559 559 560 560 ! Contribution to heat flux to the ocean [W.m-2], <0 561 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice561 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_rdtice 562 562 563 563 ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok) 564 sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdeltah(ji,jk) * rhoic* r1_rdtice564 sfx_bom_1d(ji) = sfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice 565 565 566 566 ! Total heat flux used in this process [W.m-2], >0 567 hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice567 hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_rdtice 568 568 569 569 ! Contribution to mass flux 570 wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice570 wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 571 571 572 572 ! update heat content (J.m-2) and layer thickness … … 576 576 577 577 ENDIF 578 END DO ! ji 579 END DO ! jk 580 581 !------------------------------------------------------------------------------! 582 ! Excessive ablation in a 1-category model 583 ! in a 1-category sea ice model, bottom ablation must not exceed hmelt (-0.15) 584 !------------------------------------------------------------------------------! 585 ! ??? keep ??? 586 ! clem bug: I think this should be included above, so we would not have to 587 ! track heat/salt/mass fluxes backwards 588 ! IF( jpl == 1 ) THEN 589 ! DO ji = kideb, kiut 590 ! IF( zf_tt(ji) >= 0._wp ) THEN 591 ! zdh = MAX( hmelt , dh_i_bott(ji) ) 592 ! zdvres = zdh - dh_i_bott(ji) ! >=0 593 ! dh_i_bott(ji) = zdh 594 ! 595 ! ! excessive energy is sent to lateral ablation 596 ! rswitch = MAX( 0._wp, SIGN( 1._wp , 1._wp - at_i_1d(ji) - epsi20 ) ) 597 ! zq_1cat(ji) = rswitch * rhoic * lfus * at_i_1d(ji) / MAX( 1._wp - at_i_1d(ji) , epsi20 ) * zdvres ! J.m-2 >=0 598 ! 599 ! ! correct salt and mass fluxes 600 ! sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_1d(ji) * a_i_1d(ji) * zdvres * rhoic * r1_rdtice ! this is only a raw approximation 601 ! wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_1d(ji) * zdvres * r1_rdtice 602 ! ENDIF 603 ! END DO 604 ! ENDIF 578 END DO 579 END DO 605 580 606 581 !------------------------------------------- … … 619 594 DO ji = kideb, kiut 620 595 zq_rema(ji) = zq_su(ji) + zq_bo(ji) 621 ! zindh = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) ) ! =1 if snow 622 ! zindq = 1._wp - MAX( 0._wp, SIGN( 1._wp, - zq_s(ji) + epsi20 ) ) 623 ! zdeltah (ji,1) = - zindh * zindq * zq_rema(ji) / MAX( zq_s(ji), epsi20 ) 624 ! zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_1d(ji) ) ) ! bound melting 625 ! zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,1) 626 ! dh_s_tot (ji) = dh_s_tot(ji) + zdeltah(ji,1) 627 ! ht_s_1d (ji) = ht_s_1d(ji) + zdeltah(ji,1) 628 ! 629 ! zq_rema(ji) = zq_rema(ji) + zdeltah(ji,1) * zq_s(ji) ! update available heat (J.m-2) 630 ! ! heat used to melt snow 631 ! hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zq_s(ji) * r1_rdtice ! W.m-2 (>0) 632 ! ! Contribution to mass flux 633 ! wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice 634 ! 596 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_1d(ji) ) ) ! =1 if snow 597 rswitch = rswitch * MAX( 0._wp, SIGN( 1._wp, q_s_1d(ji,1) - epsi20 ) ) 598 zdeltah (ji,1) = - rswitch * zq_rema(ji) / MAX( q_s_1d(ji,1), epsi20 ) 599 zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_1d(ji) ) ) ! bound melting 600 dh_s_tot (ji) = dh_s_tot(ji) + zdeltah(ji,1) 601 ht_s_1d (ji) = ht_s_1d(ji) + zdeltah(ji,1) 602 603 zq_rema(ji) = zq_rema(ji) + zdeltah(ji,1) * q_s_1d(ji,1) ! update available heat (J.m-2) 604 ! heat used to melt snow 605 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * q_s_1d(ji,1) * r1_rdtice ! W.m-2 (>0) 606 ! Contribution to mass flux 607 wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice 608 ! 635 609 ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 636 610 ! Remaining heat flux (W.m-2) is sent to the ocean heat budget 637 hfx_out(ii,ij) = hfx_out(ii,ij) + ( zq_ 1cat(ji) + zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice638 639 IF( ln_ nicep.AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji)611 hfx_out(ii,ij) = hfx_out(ii,ij) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_rdtice 612 613 IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji) 640 614 END DO 641 615 … … 650 624 dh_snowice(ji) = MAX( 0._wp , ( rhosn * ht_s_1d(ji) + (rhoic-rau0) * ht_i_1d(ji) ) / ( rhosn+rau0-rhoic ) ) 651 625 652 ht_i_1d(ji) 653 ht_s_1d(ji) 626 ht_i_1d(ji) = ht_i_1d(ji) + dh_snowice(ji) 627 ht_s_1d(ji) = ht_s_1d(ji) - dh_snowice(ji) 654 628 655 629 ! Salinity of snow ice 656 630 ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 657 zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) /rhoic + ( 1. - zswitch_sal ) * sm_i_1d(ji)631 zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) * r1_rhoic + ( 1. - zswitch_sal ) * sm_i_1d(ji) 658 632 659 633 ! entrapment during snow ice formation 660 ! new salinity difference stored (to be used in limthd_ ent.F90)661 IF ( n um_sal == 2 ) THEN662 rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi 10 ) )634 ! new salinity difference stored (to be used in limthd_sal.F90) 635 IF ( nn_icesal == 2 ) THEN 636 rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) ) 663 637 ! salinity dif due to snow-ice formation 664 dsm_i_si_1d(ji) = ( zs_snic - sm_i_1d(ji) ) * dh_snowice(ji) / MAX( ht_i_1d(ji), epsi 10 ) * rswitch638 dsm_i_si_1d(ji) = ( zs_snic - sm_i_1d(ji) ) * dh_snowice(ji) / MAX( ht_i_1d(ji), epsi20 ) * rswitch 665 639 ! salinity dif due to bottom growth 666 640 IF ( zf_tt(ji) < 0._wp ) THEN 667 dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_1d(ji) ) * dh_i_bott(ji) / MAX( ht_i_1d(ji), epsi 10 ) * rswitch641 dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_1d(ji) ) * dh_i_bott(ji) / MAX( ht_i_1d(ji), epsi20 ) * rswitch 668 642 ENDIF 669 643 ENDIF … … 691 665 h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji) 692 666 693 ! Total ablation (to debug) 694 IF( ht_i_1d(ji) <= 0._wp ) a_i_1d(ji) = 0._wp 695 696 END DO !ji 667 END DO 697 668 698 669 ! … … 700 671 ! Update temperature, energy 701 672 !------------------------------------------- 702 !clem bug: we should take snow into account here703 673 DO ji = kideb, kiut 704 674 rswitch = 1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_1d(ji) ) ) 705 t_su_1d(ji) = rswitch * t_su_1d(ji) + ( 1.0 - rswitch ) * rt t706 END DO ! ji675 t_su_1d(ji) = rswitch * t_su_1d(ji) + ( 1.0 - rswitch ) * rt0 676 END DO 707 677 708 678 DO jk = 1, nlay_s 709 679 DO ji = kideb,kiut 710 680 ! mask enthalpy 711 rswitch = MAX( 0._wp , SIGN( 1._wp, - ht_s_1d(ji) ) )712 q_s_1d(ji,jk) = ( 1.0 - rswitch )* q_s_1d(ji,jk)681 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp, - ht_s_1d(ji) ) ) 682 q_s_1d(ji,jk) = rswitch * q_s_1d(ji,jk) 713 683 ! recalculate t_s_1d from q_s_1d 714 t_s_1d(ji,jk) = rt t + ( 1._wp - rswitch )* ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic )684 t_s_1d(ji,jk) = rt0 + rswitch * ( - q_s_1d(ji,jk) / ( rhosn * cpic ) + lfus / cpic ) 715 685 END DO 716 686 END DO 717 687 718 CALL wrk_dealloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_1cat, zq_rema ) 719 CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s ) 720 CALL wrk_dealloc( jpij, nlay_i+1, zdeltah, zh_i ) 721 CALL wrk_dealloc( jpij, icount ) 688 ! --- ensure that a_i = 0 where ht_i = 0 --- 689 WHERE( ht_i_1d == 0._wp ) a_i_1d = 0._wp 690 691 CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema ) 692 CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s, zsnw ) 693 CALL wrk_dealloc( jpij, nlay_i, zdeltah, zh_i ) 694 CALL wrk_dealloc( jpij, nlay_i, icount ) 722 695 ! 723 696 ! 724 697 END SUBROUTINE lim_thd_dh 698 699 700 !!-------------------------------------------------------------------------- 701 !! INTERFACE lim_thd_snwblow 702 !! ** Purpose : Compute distribution of precip over the ice 703 !!-------------------------------------------------------------------------- 704 SUBROUTINE lim_thd_snwblow_2d( pin, pout ) 705 REAL(wp), DIMENSION(:,:), INTENT(in) :: pin ! previous fraction lead ( pfrld or (1. - a_i_b) ) 706 REAL(wp), DIMENSION(:,:), INTENT(out) :: pout 707 pout = ( 1._wp - ( pin )**rn_betas ) 708 END SUBROUTINE lim_thd_snwblow_2d 709 710 SUBROUTINE lim_thd_snwblow_1d( pin, pout ) 711 REAL(wp), DIMENSION(:), INTENT(in) :: pin 712 REAL(wp), DIMENSION(:), INTENT(out) :: pout 713 pout = ( 1._wp - ( pin )**rn_betas ) 714 END SUBROUTINE lim_thd_snwblow_1d 715 725 716 726 717 #else
Note: See TracChangeset
for help on using the changeset viewer.
