7 !! Computation of the surface temperature (must be less than 0 deg C!)
8 !! and of the accumulation-ablation function.
12 !! Copyright 2009-2013 Ralf Greve
16 !! This file is part of SICOPOLIS.
18 !! SICOPOLIS is free software: you can redistribute it and/or modify
19 !! it under the terms of the GNU General Public License as published by
20 !! the Free Software Foundation, either version 3 of the License, or
21 !! (at your option) any later version.
23 !! SICOPOLIS is distributed in the hope that it will be useful,
24 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
25 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
26 !! GNU General Public License for more details.
28 !! You should have received a copy of the GNU General Public License
29 !! along with SICOPOLIS. If not, see <http://www.gnu.org/licenses/>.
31 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
34 !> Computation of the surface temperature (must be less than 0 deg C!)
35 !! and of the accumulation-ablation function.
36 !<------------------------------------------------------------------------------
38 delta_ts, glac_index, z_sl, dzsl_dtau, z_mar)
48 real(dp),
intent(in) :: time, dtime, dxi, deta
50 real(dp),
intent(out) :: delta_ts, glac_index, dzsl_dtau, z_mar
51 real(dp),
intent(inout) :: z_sl
59 integer(i4b) :: i, j, n
60 integer(i4b) :: i_gr, i_kl
61 integer(i4b) :: nrec_temp_precip
62 integer(i4b),
save :: nrec_temp_precip_save = -1
64 real(dp) :: z_sl_min, t1, t2, t3, t4, t5, t6
65 real(dp) :: time_gr, time_kl
66 real(dp) :: z_sle_present, z_sle_help
67 real(dp),
dimension(0:JMAX,0:IMAX,0:12) :: precip
68 real(dp),
dimension(0:JMAX,0:IMAX) :: &
69 snowfall, rainfall, melt, melt_star
70 real(dp),
dimension(0:JMAX,0:IMAX,12) :: temp_mm
71 real(dp),
dimension(0:JMAX,0:IMAX) :: temp_ma, temp_ampl
72 real(dp),
dimension(0:JMAX,0:IMAX) :: temp_ma_anom, temp_mj_anom, precip_ma_anom
73 real(dp),
dimension(0:IMAX,0:JMAX),
save :: temp_ma_anom_tra, temp_mj_anom_tra, &
78 real(dp),
dimension(12) :: temp_mm_help
79 real(dp) :: temp_jja_help
80 real(dp),
dimension(0:JMAX,0:IMAX) :: et
81 real(dp) :: theta_ma, c_ma, kappa_ma, gamma_ma, &
82 theta_mj, c_mj, kappa_mj, gamma_mj
83 real(dp) :: sine_factor
84 real(dp) :: gamma_p, zs_thresh, &
85 temp_rain, temp_snow, &
86 inv_delta_temp_rain_snow, coeff(0:5), inv_sqrt2_s_stat, &
87 precip_fact, frac_solid
89 phi_sep, temp_lt, temp_ht, inv_delta_temp_ht_lt, &
90 beta1_lt, beta1_ht, beta2_lt, beta2_ht, &
91 beta1, beta2, pmax, mu, lambda_lti, temp_lti
92 real(dp),
dimension(256) :: pdd_mod_lat, delta_pdd_mod_lat_inv, &
93 pdd_mod_fac_w, pdd_mod_fac_e, pdd_mod_fac
94 real(dp) :: lon_w_e_sep, pdd_mod_fac_interpol
96 logical,
dimension(0:JMAX,0:IMAX) :: check_point
101 integer(i4b) :: nc3cor(3)
103 integer(i4b) :: nc3cnt(3)
107 real(dp),
parameter :: &
108 inv_twelve = 1.0_dp/12.0_dp, one_third = 1.0_dp/3.0_dp
127 delta_ts = sine_amplit &
128 *cos(2.0_dp*pi*time/(sine_period*year_sec)) &
135 if (time/year_sec <
real(grip_time_min,dp)) then
136 delta_ts = griptemp(0)
137 else if (time/year_sec <
real(grip_time_max,dp)) then
139 i_kl = floor(((time/year_sec) &
140 -
real(grip_time_min,dp))/
real(grip_time_stp,dp))
143 i_gr = ceiling(((time/year_sec) &
144 -
real(grip_time_min,dp))/
real(grip_time_stp,dp))
145 i_gr = min(i_gr, ndata_grip)
147 if (i_kl == i_gr)
then
149 delta_ts = griptemp(i_kl)
153 time_kl = (grip_time_min + i_kl*grip_time_stp) *year_sec
154 time_gr = (grip_time_min + i_gr*grip_time_stp) *year_sec
156 delta_ts = griptemp(i_kl) &
157 +(griptemp(i_gr)-griptemp(i_kl)) &
158 *(time-time_kl)/(time_gr-time_kl)
164 delta_ts = griptemp(ndata_grip)
167 delta_ts = delta_ts * grip_temp_fact
174 if (time/year_sec <
real(gi_time_min,dp)) then
175 glac_index = glacial_index(0)
176 else if (time/year_sec <
real(gi_time_max,dp)) then
178 i_kl = floor(((time/year_sec) &
179 -
real(gi_time_min,dp))/
real(gi_time_stp,dp))
182 i_gr = ceiling(((time/year_sec) &
183 -
real(gi_time_min,dp))/
real(gi_time_stp,dp))
184 i_gr = min(i_gr, ndata_gi)
186 if (i_kl == i_gr)
then
188 glac_index = glacial_index(i_kl)
192 time_kl = (gi_time_min + i_kl*gi_time_stp) *year_sec
193 time_gr = (gi_time_min + i_gr*gi_time_stp) *year_sec
195 glac_index = glacial_index(i_kl) &
196 +(glacial_index(i_gr)-glacial_index(i_kl)) &
197 *(time-time_kl)/(time_gr-time_kl)
203 glac_index = glacial_index(ndata_gi)
209 #elif ( (TSURFACE==6) && (ACCSURFACE==6) )
211 if (time/year_sec <= temp_precip_time_min)
then
213 else if (time/year_sec < temp_precip_time_max)
then
214 nrec_temp_precip = nint( ((time/year_sec)-temp_precip_time_min) &
215 / temp_precip_time_stp )
217 nrec_temp_precip = ndata_temp_precip
220 if ( nrec_temp_precip < 0 )
then
221 stop
' boundary: nrec_temp_precip < 0, not allowed!'
222 else if ( nrec_temp_precip > ndata_temp_precip )
then
223 stop
' boundary: nrec_temp_precip > ndata_temp_precip, not allowed!'
226 if ( nrec_temp_precip /= nrec_temp_precip_save )
then
230 nc3cor(3) = nrec_temp_precip + 1
235 call
check( nf90_inq_varid(ncid_temp_precip,
'annualtemp_anom', ncv) )
236 call
check( nf90_get_var(ncid_temp_precip, ncv, temp_ma_anom_tra, &
237 start=nc3cor, count=nc3cnt) )
239 call
check( nf90_inq_varid(ncid_temp_precip,
'julytemp_anom', ncv) )
240 call
check( nf90_get_var(ncid_temp_precip, ncv, temp_mj_anom_tra, &
241 start=nc3cor, count=nc3cnt) )
243 call
check( nf90_inq_varid(ncid_temp_precip,
'precipitation_anom', ncv) )
244 call
check( nf90_get_var(ncid_temp_precip, ncv, precip_ma_anom_tra, &
245 start=nc3cor, count=nc3cnt) )
249 temp_ma_anom = transpose(temp_ma_anom_tra)
250 temp_mj_anom = transpose(temp_mj_anom_tra)
251 precip_ma_anom = transpose(precip_ma_anom_tra) *(1.0e-03_dp/year_sec)*(rho_w/rho)
254 nrec_temp_precip_save = nrec_temp_precip
269 t1 = -250000.0_dp *year_sec
270 t2 = -140000.0_dp *year_sec
271 t3 = -125000.0_dp *year_sec
272 t4 = -21000.0_dp *year_sec
273 t5 = -8000.0_dp *year_sec
274 t6 = 0.0_dp *year_sec
278 else if (time < t2)
then
279 z_sl = z_sl_min*(time-t1)/(t2-t1)
280 else if (time < t3)
then
281 z_sl = -z_sl_min*(time-t3)/(t3-t2)
282 else if (time < t4)
then
283 z_sl = z_sl_min*(time-t3)/(t4-t3)
284 else if (time < t5)
then
285 z_sl = -z_sl_min*(time-t5)/(t5-t4)
286 else if (time < t6)
then
296 if (time/year_sec <
real(specmap_time_min,dp)) then
297 z_sl = specmap_zsl(0)
298 else if (time/year_sec <
real(specmap_time_max,dp)) then
300 i_kl = floor(((time/year_sec) &
301 -
real(specmap_time_min,dp))/
real(specmap_time_stp,dp))
304 i_gr = ceiling(((time/year_sec) &
305 -
real(specmap_time_min,dp))/
real(specmap_time_stp,dp))
306 i_gr = min(i_gr, ndata_specmap)
308 if (i_kl == i_gr)
then
310 z_sl = specmap_zsl(i_kl)
314 time_kl = (specmap_time_min + i_kl*specmap_time_stp) *year_sec
315 time_gr = (specmap_time_min + i_gr*specmap_time_stp) *year_sec
317 z_sl = specmap_zsl(i_kl) &
318 +(specmap_zsl(i_gr)-specmap_zsl(i_kl)) &
319 *(time-time_kl)/(time_gr-time_kl)
325 z_sl = specmap_zsl(ndata_specmap)
332 if ( z_sl_old > -999999.9_dp )
then
333 dzsl_dtau = (z_sl-z_sl_old)/dtime
342 #if ( MARINE_ICE_CALVING==2 || MARINE_ICE_CALVING==3 )
344 #elif ( MARINE_ICE_CALVING==4 || MARINE_ICE_CALVING==5 )
345 z_mar = fact_z_mar*z_sl
346 #elif ( MARINE_ICE_CALVING==6 || MARINE_ICE_CALVING==7 )
347 if (z_sl >= -80.0_dp)
then
350 z_mar = 10.25_dp*(z_sl+80.0_dp)-200.0_dp
352 z_mar = fact_z_mar*z_mar
364 check_point(j,i) = .false.
370 if (maske(j,i) >= 2)
then
371 check_point(j ,i ) = .true.
372 check_point(j ,i+1) = .true.
373 check_point(j ,i-1) = .true.
374 check_point(j+1,i ) = .true.
375 check_point(j-1,i ) = .true.
382 if (check_point(j,i))
then
392 if (check_point(j,i))
then
393 maske(j,i) = maske_neu(j,i)
400 #if TEMP_PRESENT_PARA == 1 /* Parameterization by Ritz et al. (1997) */
403 gamma_ma = -7.992e-03_dp
408 gamma_mj = -6.277e-03_dp
412 #elif TEMP_PRESENT_PARA == 2 /* Parameterization by Fausto et al. (2009) */
415 gamma_ma = -6.309e-03_dp
417 kappa_ma = -0.0672_dp
420 gamma_mj = -5.426e-03_dp
422 kappa_mj = -0.0518_dp
426 stop
' boundary: Parameter TEMP_PRESENT_PARA must be either 1 or 2!'
435 temp_ma_present(j,i) = theta_ma &
437 + c_ma*phi(j,i)*pi_180_inv &
438 + kappa_ma*(modulo(lambda(j,i)+pi,2.0_dp*pi)-pi)*pi_180_inv
443 temp_mj_present(j,i) = theta_mj &
445 + c_mj*phi(j,i)*pi_180_inv &
446 + kappa_mj*(modulo(lambda(j,i)+pi,2.0_dp*pi)-pi)*pi_180_inv
453 temp_ma(j,i) = temp_ma_present(j,i) + delta_ts
454 temp_mm(j,i,7) = temp_mj_present(j,i) + delta_ts
456 #elif (TSURFACE == 5)
460 temp_ma(j,i) = temp_ma_present(j,i) + glac_index*temp_ma_lgm_anom(j,i)
461 temp_mm(j,i,7) = temp_mj_present(j,i) + glac_index*temp_mj_lgm_anom(j,i)
463 #elif (TSURFACE == 6)
468 temp_ma(j,i) = temp_ma_present(j,i) + temp_ma_anom_fact*temp_ma_anom(j,i)
469 temp_mm(j,i,7) = temp_mj_present(j,i) + temp_mj_anom_fact*temp_mj_anom(j,i)
475 temp_ampl(j,i) = temp_mm(j,i,7) - temp_ma(j,i)
477 if (temp_ampl(j,i) < eps)
then
488 sine_factor = sin((
real(n,dp)-4.0_dp)*pi/6.0_dp)
492 temp_mm(j,i,n) = temp_ma(j,i) + sine_factor*temp_ampl(j,i)
500 #if (ELEV_DESERT == 1)
502 gamma_p = gamma_p*1.0e-03_dp
504 zs_thresh = zs_thresh
508 #if (SOLID_PRECIP == 1) /* Marsiat (1994) */
515 inv_delta_temp_rain_snow = 1.0_dp/(temp_rain-temp_snow)
517 #elif (SOLID_PRECIP == 2) /* Bales et al. (2009) */
524 coeff(0) = 5.4714e-01_dp
525 coeff(1) = -9.1603e-02_dp
526 coeff(2) = -3.314e-03_dp
527 coeff(3) = 4.66e-04_dp
528 coeff(4) = 3.8e-05_dp
529 coeff(5) = 6.0e-07_dp
531 #elif (SOLID_PRECIP == 3) /* Huybrechts and de Wolde (1999) */
535 temp_snow = temp_rain
541 inv_sqrt2_s_stat = 1.0_dp/(sqrt(2.0_dp)*s_stat)
545 #if (ABLSURFACE==1 || ABLSURFACE==2)
549 phi_sep = phi_sep_0*pi_180
554 inv_delta_temp_ht_lt = 1.0_dp/(temp_ht-temp_lt)
556 beta1_lt = beta1_lt_0 *(0.001_dp/86400.0_dp)*(rho_w/rho)
558 beta1_ht = beta1_ht_0 *(0.001_dp/86400.0_dp)*(rho_w/rho)
560 beta2_lt = beta2_lt_0 *(0.001_dp/86400.0_dp)*(rho_w/rho)
562 beta2_ht = beta2_ht_0 *(0.001_dp/86400.0_dp)*(rho_w/rho)
565 mu = mu_0 *(1000.0_dp*86400.0_dp)*(rho/rho_w)
567 #if (PDD_MODIFIER==2)
569 lon_w_e_sep = lon_w_e_sep *pi_180
571 if (lon_w_e_sep < 0.0_dp)
then
572 lon_w_e_sep = lon_w_e_sep + 2.0_dp*pi
573 else if (lon_w_e_sep >= (2.0_dp*pi))
then
574 lon_w_e_sep = lon_w_e_sep - 2.0_dp*pi
578 pdd_mod_fac_w = 0.0_dp
579 pdd_mod_fac_e = 0.0_dp
581 pdd_mod_lat( 1) = pdd_mod_lat_01 *pi_180
582 pdd_mod_fac_w( 1) = pdd_mod_fac_w_01
583 pdd_mod_fac_e( 1) = pdd_mod_fac_e_01
584 pdd_mod_lat( 2) = pdd_mod_lat_02 *pi_180
585 pdd_mod_fac_w( 2) = pdd_mod_fac_w_02
586 pdd_mod_fac_e( 2) = pdd_mod_fac_e_02
587 pdd_mod_lat( 3) = pdd_mod_lat_03 *pi_180
588 pdd_mod_fac_w( 3) = pdd_mod_fac_w_03
589 pdd_mod_fac_e( 3) = pdd_mod_fac_e_03
590 pdd_mod_lat( 4) = pdd_mod_lat_04 *pi_180
591 pdd_mod_fac_w( 4) = pdd_mod_fac_w_04
592 pdd_mod_fac_e( 4) = pdd_mod_fac_e_04
593 pdd_mod_lat( 5) = pdd_mod_lat_05 *pi_180
594 pdd_mod_fac_w( 5) = pdd_mod_fac_w_05
595 pdd_mod_fac_e( 5) = pdd_mod_fac_e_05
596 pdd_mod_lat( 6) = pdd_mod_lat_06 *pi_180
597 pdd_mod_fac_w( 6) = pdd_mod_fac_w_06
598 pdd_mod_fac_e( 6) = pdd_mod_fac_e_06
599 pdd_mod_lat( 7) = pdd_mod_lat_07 *pi_180
600 pdd_mod_fac_w( 7) = pdd_mod_fac_w_07
601 pdd_mod_fac_e( 7) = pdd_mod_fac_e_07
602 pdd_mod_lat( 8) = pdd_mod_lat_08 *pi_180
603 pdd_mod_fac_w( 8) = pdd_mod_fac_w_08
604 pdd_mod_fac_e( 8) = pdd_mod_fac_e_08
605 pdd_mod_lat( 9) = pdd_mod_lat_09 *pi_180
606 pdd_mod_fac_w( 9) = pdd_mod_fac_w_09
607 pdd_mod_fac_e( 9) = pdd_mod_fac_e_09
608 pdd_mod_lat( 10) = pdd_mod_lat_10 *pi_180
609 pdd_mod_fac_w(10) = pdd_mod_fac_w_10
610 pdd_mod_fac_e(10) = pdd_mod_fac_e_10
612 delta_pdd_mod_lat_inv = 0.0_dp
615 delta_pdd_mod_lat_inv(n) = 1.0_dp/(pdd_mod_lat(n+1)-pdd_mod_lat(n))
620 #elif (ABLSURFACE==3)
622 lambda_lti = lambda_lti *(0.001_dp/year_sec)*(rho_w/rho)
632 #if (ABLSURFACE==1 || ABLSURFACE==2)
634 if (phi(j,i) <= phi_sep)
then
641 if (temp_mm(j,i,7) >= temp_ht)
then
644 else if (temp_mm(j,i,7) <= temp_lt)
then
649 + (beta1_ht-beta1_lt) &
650 *inv_delta_temp_ht_lt*(temp_mm(j,i,7)-temp_lt)
652 + (beta2_lt-beta2_ht) &
653 *(inv_delta_temp_ht_lt*(temp_ht-temp_mm(j,i,7)))**3
658 #if (PDD_MODIFIER==2)
660 if ( lambda(j,i) <= lon_w_e_sep )
then
661 pdd_mod_fac = pdd_mod_fac_w
663 pdd_mod_fac = pdd_mod_fac_e
666 if (phi(j,i) <= pdd_mod_lat(1))
then
668 beta1 = beta1 * pdd_mod_fac(1)
669 beta2 = beta2 * pdd_mod_fac(1)
671 else if (phi(j,i) >= pdd_mod_lat(n_pdd_mod))
then
673 beta1 = beta1 * pdd_mod_fac(n_pdd_mod)
674 beta2 = beta2 * pdd_mod_fac(n_pdd_mod)
680 if ( (phi(j,i) >= pdd_mod_lat(n)) &
682 (phi(j,i) <= pdd_mod_lat(n+1)) )
then
684 pdd_mod_fac_interpol = pdd_mod_fac(n) &
685 + (pdd_mod_fac(n+1)-pdd_mod_fac(n)) &
686 *(phi(j,i)-pdd_mod_lat(n)) &
687 *delta_pdd_mod_lat_inv(n)
689 beta1 = beta1 * pdd_mod_fac_interpol
690 beta2 = beta2 * pdd_mod_fac_interpol
706 #if (ACCSURFACE <= 3)
710 #if (ELEV_DESERT == 0)
714 #elif (ELEV_DESERT == 1)
716 if (zs_ref(j,i) < zs_thresh)
then
718 = exp(gamma_p*(max(zs(j,i),zs_thresh)-zs_thresh))
721 = exp(gamma_p*(max(zs(j,i),zs_thresh)-zs_ref(j,i)))
725 stop
' boundary: Parameter ELEV_DESERT must be either 0 or 1!'
729 precip(j,i,n) = precip_present(j,i,n)*precip_fact
737 precip_fact = accfact
739 precip_fact = 1.0_dp + gamma_s*delta_ts
741 precip_fact = exp(gamma_s*delta_ts)
744 #if (ACCSURFACE <= 3)
746 precip(j,i,0) = 0.0_dp
749 precip(j,i,n) = precip(j,i,n)*precip_fact
750 precip(j,i,0) = precip(j,i,0) + precip(j,i,n)*inv_twelve
754 #elif (ACCSURFACE == 5)
756 precip(j,i,0) = 0.0_dp
760 #if (PRECIP_ANOM_INTERPOL==1)
761 precip_fact = 1.0_dp-glac_index+glac_index*precip_lgm_anom(j,i,n)
763 #elif (PRECIP_ANOM_INTERPOL==2)
764 precip_fact = exp(-glac_index*gamma_precip_lgm_anom(j,i,n))
768 precip(j,i,n) = precip_present(j,i,n)*precip_fact
769 precip(j,i,0) = precip(j,i,0) + precip(j,i,n)*inv_twelve
773 #elif (ACCSURFACE == 6)
777 precip(j,i,0) = precip_ma_present(j,i) + precip_anom_fact*precip_ma_anom(j,i)
781 precip(j,i,n) = precip(j,i,0)
789 accum(j,i) = precip(j,i,0)
791 snowfall(j,i) = 0.0_dp
795 #if (SOLID_PRECIP == 1) /* Marsiat (1994) */
797 if (temp_mm(j,i,n) >= temp_rain)
then
799 else if (temp_mm(j,i,n) <= temp_snow)
then
802 frac_solid = (temp_rain-temp_mm(j,i,n))*inv_delta_temp_rain_snow
805 #elif (SOLID_PRECIP == 2) /* Bales et al. (2009) */
807 if (temp_mm(j,i,n) >= temp_rain)
then
809 else if (temp_mm(j,i,n) <= temp_snow)
then
812 frac_solid = coeff(0) + temp_mm(j,i,n) * ( coeff(1) &
813 + temp_mm(j,i,n) * ( coeff(2) &
814 + temp_mm(j,i,n) * ( coeff(3) &
815 + temp_mm(j,i,n) * ( coeff(4) &
816 + temp_mm(j,i,n) * coeff(5) ) ) ) )
820 #elif (SOLID_PRECIP == 3) /* Huybrechts and de Wolde (1999) */
822 frac_solid = 1.0_dp &
823 - 0.5_dp*
erfcc((temp_rain-temp_mm(j,i,n))*inv_sqrt2_s_stat)
827 snowfall(j,i) = snowfall(j,i) + precip(j,i,n)*frac_solid*inv_twelve
831 rainfall(j,i) = precip(j,i,0) - snowfall(j,i)
833 if (snowfall(j,i) < 0.0_dp) snowfall(j,i) = 0.0_dp
834 if (rainfall(j,i) < 0.0_dp) rainfall(j,i) = 0.0_dp
840 #if (ABLSURFACE==1 || ABLSURFACE==2)
845 temp_mm_help(n) = temp_mm(j,i,n)
848 call
pdd(temp_mm_help, s_stat, et(j,i))
855 if ((beta1*et(j,i)) <= (pmax*snowfall(j,i)))
then
856 melt_star(j,i) = beta1*et(j,i)
858 runoff(j,i) = melt(j,i)+rainfall(j,i)
860 melt_star(j,i) = pmax*snowfall(j,i)
861 melt(j,i) = beta2*(et(j,i)-melt_star(j,i)/beta1)
862 runoff(j,i) = melt(j,i)+rainfall(j,i)
865 #elif (ABLSURFACE==2)
867 if ( rainfall(j,i) <= (pmax*snowfall(j,i)) )
then
869 if ( (rainfall(j,i)+beta1*et(j,i)) <= (pmax*snowfall(j,i)) )
then
870 melt_star(j,i) = rainfall(j,i)+beta1*et(j,i)
872 runoff(j,i) = melt(j,i)
874 melt_star(j,i) = pmax*snowfall(j,i)
876 *(et(j,i)-(melt_star(j,i)-rainfall(j,i))/beta1)
877 runoff(j,i) = melt(j,i)
882 melt_star(j,i) = pmax*snowfall(j,i)
883 melt(j,i) = beta2*et(j,i)
884 runoff(j,i) = melt(j,i) + rainfall(j,i)-pmax*snowfall(j,i)
890 #elif (ABLSURFACE==3)
892 temp_jja_help = one_third*(temp_mm(j,i,6)+temp_mm(j,i,7)+temp_mm(j,i,8))
894 melt_star(j,i) = 0.0_dp
895 melt(j,i) = lambda_lti*max((temp_jja_help-temp_lti), 0.0_dp)
896 runoff(j,i) = melt(j,i) + rainfall(j,i)
906 as_perp(j,i) = accum(j,i) - evap(j,i) - runoff(j,i)
912 if (melt_star(j,i) >= melt(j,i))
then
913 temp_s(j,i) = temp_ma(j,i) &
914 +mu*(melt_star(j,i)-melt(j,i))
916 temp_s(j,i) = temp_ma(j,i)
919 if (temp_s(j,i) > -0.001_dp) temp_s(j,i) = -0.001_dp