44 #if !defined(ALLOW_OPENAD) /* Normal */ 62 real(dp),
intent(in) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
63 real(dp),
intent(in) :: dtime_temp
65 integer(i4b) :: i, j, kc, kt, kr, ii, jj
66 real(dp) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
67 at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
68 at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
69 acb1, acb2, acb3, acb4, &
70 ai1(0:kcmax), ai2(0:kcmax), ai3, &
72 real(dp) :: aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld
73 real(dp) :: dtime_temp_inv, dtt_2dxi, dtt_2deta
74 real(dp) :: temp_c_help(0:kcmax)
75 real(dp) :: fact_thick
76 real(dp) :: time_lag_cts
77 real(dp) :: Vol_t, Vol_t_smooth, korrfakt_t, &
78 dH_t_smooth(0:jmax,0:imax)
82 at7 = 2.0_dp/
rho*dtime_temp
84 aw1 = dtime_temp/dzeta_t
85 aw2 = dtime_temp/dzeta_t
86 aw3 = dtime_temp/dzeta_t
87 aw4 = dtime_temp/dzeta_t
88 aw5 =
nue/
rho*dtime_temp/(dzeta_t**2)
89 aw7 = 2.0_dp/(
rho*
l)*dtime_temp
92 aw9 =
beta/
l*dtime_temp
94 ai3 = agediff*dtime_temp/(dzeta_t**2)
107 acb1 = (
ea-1.0_dp)/
aa/dzeta_c
109 acb1 = 1.0_dp/dzeta_c
118 aqtld = dzeta_t/dtime_temp
120 dtt_2dxi = 0.5_dp*dtime_temp/dxi
121 dtt_2deta = 0.5_dp*dtime_temp/deta
123 dtime_temp_inv = 1.0_dp/dtime_temp
129 at1(kc) = (
ea-1.0_dp)/(
aa*
eaz_c(kc))*dtime_temp/dzeta_c
130 at2_1(kc) = (
ea-1.0_dp)/(
aa*
eaz_c(kc))*dtime_temp/dzeta_c
133 at3_1(kc) = (
ea-1.0_dp)/(
aa*
eaz_c(kc))*dtime_temp/dzeta_c
136 at4_1(kc) = (
ea-1.0_dp)/(
aa*
eaz_c(kc))*dtime_temp/dzeta_c
141 if (kc /= kcmax)
then 142 at6(kc) = (
ea-1.0_dp) &
148 ai1(kc) = agediff*(
ea-1.0_dp)/(
aa*
eaz_c(kc)) &
150 if (kc /= kcmax)
then 151 ai2(kc) = (
ea-1.0_dp) &
160 at1(kc) = dtime_temp/dzeta_c
161 at2_1(kc) = dtime_temp/dzeta_c
164 at3_1(kc) = dtime_temp/dzeta_c
167 at4_1(kc) = dtime_temp/dzeta_c
170 at5(kc) = 1.0_dp/
rho &
172 if (kc /= kcmax)
then 180 if (kc /= kcmax)
then 196 if (
maske(j,i)==0_i1b)
then 200 if (
n_cts(j,i) == -1_i1b)
then 207 call calc_temp1(at1, at2_1, at2_2, at3_1, at3_2, &
208 at4_1, at4_2, at5, at6, at7, atr1, acb1, acb2, &
209 acb3, acb4, alb1, ai1, ai2, &
210 dtime_temp, dtt_2dxi, dtt_2deta, i, j)
218 call calc_temp2(at1, at2_1, at2_2, at3_1, at3_2, &
219 at4_1, at4_2, at5, at6, at7, atr1, alb1, &
221 dtime_temp, dtt_2dxi, dtt_2deta, i, j)
239 call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
240 at4_1, at4_2, at5, at6, at7, atr1, &
242 aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
243 ai1, ai2, ai3, dzeta_t, &
244 dtime_temp, dtt_2dxi, dtt_2deta, i, j)
247 at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, &
248 aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
249 ai1, ai2, ai3, dzeta_t, &
250 dtime_temp, dtt_2dxi, dtt_2deta, dtime_temp_inv, &
257 else if (
n_cts(j,i) == 0_i1b)
then 264 call calc_temp2(at1, at2_1, at2_2, at3_1, at3_2, &
265 at4_1, at4_2, at5, at6, at7, atr1, alb1, &
267 dtime_temp, dtt_2dxi, dtt_2deta, i, j)
275 call calc_temp1(at1, at2_1, at2_2, at3_1, at3_2, &
276 at4_1, at4_2, at5, at6, at7, atr1, acb1, acb2, &
277 acb3, acb4, alb1, ai1, ai2, &
278 dtime_temp, dtt_2dxi, dtt_2deta, i, j)
284 call calc_temp2(at1, at2_1, at2_2, at3_1, at3_2, &
285 at4_1, at4_2, at5, at6, at7, atr1, alb1, &
287 dtime_temp, dtt_2dxi, dtt_2deta, i, j)
307 call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
308 at4_1, at4_2, at5, at6, at7, atr1, &
310 aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
311 ai1, ai2, ai3, dzeta_t, &
312 dtime_temp, dtt_2dxi, dtt_2deta, i, j)
315 at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, &
316 aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
317 ai1, ai2, ai3, dzeta_t, &
318 dtime_temp, dtt_2dxi, dtt_2deta, dtime_temp_inv, &
332 call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
333 at4_1, at4_2, at5, at6, at7, atr1, &
335 aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
336 ai1, ai2, ai3, dzeta_t, &
337 dtime_temp, dtt_2dxi, dtt_2deta, i, j)
342 at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, &
343 aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
344 ai1, ai2, ai3, dzeta_t, &
345 dtime_temp, dtt_2dxi, dtt_2deta, dtime_temp_inv, &
348 call shift_cts_downward(at1, at2_1, at2_2, at3_1, at3_2, &
349 at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, &
350 aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
351 ai1, ai2, ai3, dzeta_t, &
352 dtime_temp, dtt_2dxi, dtt_2deta, dtime_temp_inv, &
360 else if (
maske(j,i)==3_i1b)
then 367 call calc_temp_ssa(at1, at2_1, at2_2, at3_1, at3_2, &
368 at4_1, at4_2, at5, at6, at7, atr1, alb1, &
370 dtime_temp, dtt_2dxi, dtt_2deta, i, j)
403 if ( (
maske(j,i) == 0_i1b).or.(
maske(j,i) == 3_i1b) )
then 443 if ( (
maske(j,i) == 0_i1b).or.(
maske(j,i) == 3_i1b) )
then 483 if ( (
maske(j,i) == 0_i1b).or.(
maske(j,i) == 3_i1b) )
then 523 if ( (
maske(j,i) == 0_i1b).or.(
maske(j,i) == 3_i1b) )
then 566 if ( (
maske(j,i) == 0_i1b).or.(
maske(j,i) == 3_i1b) )
then 605 if ( (
maske(j,i) == 0_i1b).or.(
maske(j,i) == 3_i1b) )
then 650 if ( (
maske(j,i) == 0_i1b).or.(
maske(j,i) == 3_i1b) )
then 689 if ( (
maske(j,i) == 0_i1b).or.(
maske(j,i) == 3_i1b) )
then 750 dh_t_smooth(j,i) = numdiff_h_t* ( -4.0_dp*
h_t_neu(j,i) &
753 if (dh_t_smooth(j,i) > 0.001_dp)
n_cts_neu(j,i) = 1_i1b
769 vol_t_smooth = 0.0_dp
773 vol_t_smooth = vol_t_smooth +
h_t_neu(j,i)*
area(j,i)
780 if (vol_t_smooth > 0.0_dp)
then 782 korrfakt_t = vol_t/vol_t_smooth
797 time_lag_cts = tau_cts*year_sec
806 /(time_lag_cts+dtime_temp)
828 real(dp),
intent(in) :: dxi, deta, dzeta_c, dzeta_t, dzeta_r
829 real(dp),
intent(in) :: dtime_temp
831 integer(i4b) :: i, j, kc, kr, ii, jj
832 real(dp) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
833 at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
834 at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
835 acb1, acb2, acb3, acb4, &
836 ai1(0:kcmax), ai2(0:kcmax), ai3, &
838 real(dp) :: aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld
839 real(dp) :: dtime_temp_inv, dtt_2dxi, dtt_2deta
840 real(dp) :: temp_c_help(0:kcmax)
844 at7 = 2.0_dp/
rho*dtime_temp
846 aw1 = dtime_temp/dzeta_t
847 aw2 = dtime_temp/dzeta_t
848 aw3 = dtime_temp/dzeta_t
849 aw4 = dtime_temp/dzeta_t
850 aw5 =
nue/
rho*dtime_temp/(dzeta_t**2)
851 aw7 = 2.0_dp/(
rho*
l)*dtime_temp
854 aw9 =
beta/
l*dtime_temp
856 ai3 = agediff*dtime_temp/(dzeta_t**2)
869 acb1 = (
ea-1.0_dp)/
aa/dzeta_c
871 acb1 = 1.0_dp/dzeta_c
880 aqtld = dzeta_t/dtime_temp
882 dtt_2dxi = 0.5_dp*dtime_temp/dxi
883 dtt_2deta = 0.5_dp*dtime_temp/deta
885 dtime_temp_inv = 1.0_dp/dtime_temp
891 at1(kc) = (
ea-1.0_dp)/(
aa*
eaz_c(kc))*dtime_temp/dzeta_c
892 at2_1(kc) = (
ea-1.0_dp)/(
aa*
eaz_c(kc))*dtime_temp/dzeta_c
895 at3_1(kc) = (
ea-1.0_dp)/(
aa*
eaz_c(kc))*dtime_temp/dzeta_c
898 at4_1(kc) = (
ea-1.0_dp)/(
aa*
eaz_c(kc))*dtime_temp/dzeta_c
903 if (kc /= kcmax)
then 904 at6(kc) = (
ea-1.0_dp) &
910 ai1(kc) = agediff*(
ea-1.0_dp)/(
aa*
eaz_c(kc)) &
912 if (kc /= kcmax)
then 913 ai2(kc) = (
ea-1.0_dp) &
922 at1(kc) = dtime_temp/dzeta_c
923 at2_1(kc) = dtime_temp/dzeta_c
926 at3_1(kc) = dtime_temp/dzeta_c
929 at4_1(kc) = dtime_temp/dzeta_c
932 at5(kc) = 1.0_dp/
rho &
934 if (kc /= kcmax)
then 942 if (kc /= kcmax)
then 958 if (
maske(j,i)==0_i1b)
then 965 call calc_temp1(at1, at2_1, at2_2, at3_1, at3_2, &
966 at4_1, at4_2, at5, at6, at7, atr1, acb1, acb2, &
967 acb3, acb4, alb1, ai1, ai2, &
968 dtime_temp, dtt_2dxi, dtt_2deta, i, j)
991 else if (
maske(j,i)==3_i1b)
then 999 call calc_temp_ssa(at1, at2_1, at2_2, at3_1, at3_2, &
1000 at4_1, at4_2, at5, at6, at7, atr1, alb1, &
1002 dtime_temp, dtt_2dxi, dtt_2deta, i, j)
1036 if ( (
maske(j,i) == 0_i1b).or.(
maske(j,i) == 3_i1b) )
then 1073 if ( (
maske(j,i) == 0_i1b).or.(
maske(j,i) == 3_i1b) )
then 1110 if ( (
maske(j,i) == 0_i1b).or.(
maske(j,i) == 3_i1b) )
then 1147 if ( (
maske(j,i) == 0_i1b).or.(
maske(j,i) == 3_i1b) )
then 1187 if ( (
maske(j,i) == 0_i1b).or.(
maske(j,i) == 3_i1b) )
then 1223 if ( (
maske(j,i) == 0_i1b).or.(
maske(j,i) == 3_i1b) )
then 1265 if ( (
maske(j,i) == 0_i1b).or.(
maske(j,i) == 3_i1b) )
then 1301 if ( (
maske(j,i) == 0_i1b).or.(
maske(j,i) == 3_i1b) )
then 1348 #if (defined(TEMP_CONST)) 1349 if ( temp_const > -
eps )
then 1350 errormsg =
' >>> calc_temp_const: TEMP_CONST must be negative!' 1368 #if (defined(AGE_CONST)) 1387 subroutine calc_temp1(at1, at2_1, at2_2, at3_1, at3_2, &
1388 at4_1, at4_2, at5, at6, at7, atr1, acb1, acb2, &
1389 acb3, acb4, alb1, ai1, ai2, &
1390 dtime_temp, dtt_2dxi, dtt_2deta, i, j)
1395 #if !defined(ALLOW_OPENAD) /* Normal */ 1399 #endif /* Normal vs. OpenAD */ 1403 #if !defined(ALLOW_OPENAD) /* Normal */ 1404 integer(i4b),
intent(in) :: i, j
1406 integer(i4b),
intent(inout) :: i, j
1407 #endif /* Normal vs. OpenAD */ 1408 real(dp),
intent(in) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
1409 at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
1410 at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
1411 ai1(0:kcmax), ai2(0:kcmax), &
1412 atr1, acb1, acb2, acb3, acb4, alb1
1413 real(dp),
intent(in) :: dtime_temp, dtt_2dxi, dtt_2deta
1415 integer(i4b) :: kc, kt, kr
1416 real(dp) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), ct4(0:kcmax), &
1417 ct5(0:kcmax), ct6(0:kcmax), ct7(0:kcmax), ctr1, &
1418 ccb1, ccb2, ccb3, ccb4, clb1
1419 real(dp) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), ct3_sg(0:kcmax), &
1420 ct4_sg(0:kcmax), adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
1421 real(dp) :: ci1(0:kcmax), ci2(0:kcmax)
1422 real(dp) :: temp_c_help(0:kcmax)
1423 real(dp) :: vx_c_help, vy_c_help
1424 real(dp) :: adv_vert_help
1425 real(dp) :: dtt_dxi, dtt_deta
1426 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
1427 lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
1428 lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
1429 lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
1430 lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
1431 real(dp),
parameter :: zero=0.0_dp
1435 if ((i == 0).or.(i == imax).or.(j == 0).or.(j == jmax))
then 1436 errormsg =
' >>> calc_temp1: Boundary points not allowed!' 1453 ccb3 = acb3*0.5_dp*(
vx_t(0,j,i)+
vx_t(0,j,i-1)) &
1455 ccb4 = acb4*0.5_dp*(
vy_t(0,j,i)+
vy_t(0,j-1,i)) &
1470 clb1 = alb1*
q_geo(j,i)
1475 ct1(kc) = at1(kc)/
h_c(j,i)*0.5_dp*(
vz_c(kc,j,i)+
vz_c(kc-1,j,i))
1479 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/
h_c(j,i)*
vz_c(kc,j,i)
1482 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/
h_c(j,i)*
vz_c(kc,j,i)
1485 #elif (ADV_VERT==2 || ADV_VERT==3) 1488 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/
h_c(j,i)*
vz_c(kc,j,i)
1495 ct2(kc) = ( at2_1(kc)*
dzm_dtau(j,i) &
1518 ct7(kc) = 2.0_dp*at7 &
1522 #if !defined(ALLOW_OPENAD) 1523 enh_c(kc,j,i), 0_i1b) &
1525 enh_c(kc,j,i), 0_i4b) &
1531 ci1(kc) = ai1(kc)/
h_c(j,i)
1538 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
1539 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
1540 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
1541 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
1542 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
1544 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
1545 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
1546 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
1547 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
1548 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
1550 #elif (ADV_VERT==2 || ADV_VERT==3) 1553 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
1554 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
1555 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
1556 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
1557 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
1563 temp_c_help(kc) = 0.5_dp*(
temp_c(kc,j,i)+
temp_c(kc+1,j,i))
1567 ci2(kc) = ai2(kc)/
h_c(j,i)
1571 dtt_dxi = 2.0_dp*dtt_2dxi
1572 dtt_deta = 2.0_dp*dtt_2deta
1580 lgs_a2(kr) = -1.0_dp
1588 lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
1590 lgs_b(kr) =
temp_r(kr,j,i)
1599 lgs_a2(kr) = -1.0_dp
1600 lgs_b(kr) = 2.0_dp*clb1
1608 lgs_a1(kr) = -(ccb1+ccb2)
1610 lgs_b(kr) = ccb3+ccb4
1616 lgs_a0(krmax+kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
1618 lgs_a1(krmax+kc) = 1.0_dp+ct5(kc)*(ct6(kc)+ct6(kc-1))
1619 lgs_a2(krmax+kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
1625 = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
1629 +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
1630 -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
1631 +ct5(kc)*(ct6(kc)+ct6(kc-1))
1633 = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
1638 adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
1641 = -max(adv_vert_help, 0.0_dp) &
1645 +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp) &
1646 +ct5(kc)*(ct6(kc)+ct6(kc-1))
1648 = min(adv_vert_help, 0.0_dp) &
1655 lgs_b(krmax+kc) =
temp_c(kc,j,i) + ct7(kc) &
1657 ( (
vx_c(kc,j,i)-abs(
vx_c(kc,j,i))) &
1660 +(
vx_c(kc,j,i-1)+abs(
vx_c(kc,j,i-1))) &
1664 ( (
vy_c(kc,j,i)-abs(
vy_c(kc,j,i))) &
1667 +(
vy_c(kc,j-1,i)+abs(
vy_c(kc,j-1,i))) &
1673 vx_c_help = 0.5_dp*(
vx_c(kc,j,i)+
vx_c(kc,j,i-1))
1674 vy_c_help = 0.5_dp*(
vy_c(kc,j,i)+
vy_c(kc,j-1,i))
1676 lgs_b(krmax+kc) =
temp_c(kc,j,i) + ct7(kc) &
1678 ( min(vx_c_help, 0.0_dp) &
1681 +max(vx_c_help, 0.0_dp) &
1685 ( min(vy_c_help, 0.0_dp) &
1688 +max(vy_c_help, 0.0_dp) &
1697 lgs_a0(krmax+kc) = 0.0_dp
1698 lgs_a1(krmax+kc) = 1.0_dp
1699 lgs_b(krmax+kc) =
temp_s(j,i)
1703 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax+krmax)
1729 lgs_a1(kc) = 1.0_dp - min(adv_vert_sg(kc), 0.0_dp)
1730 lgs_a2(kc) = min(adv_vert_sg(kc), 0.0_dp)
1734 lgs_b(kc) =
age_c(kc,j,i) + dtime_temp &
1736 ( (
vx_c(kc,j,i)-abs(
vx_c(kc,j,i))) &
1739 +(
vx_c(kc,j,i-1)+abs(
vx_c(kc,j,i-1))) &
1743 ( (
vy_c(kc,j,i)-abs(
vy_c(kc,j,i))) &
1746 +(
vy_c(kc,j-1,i)+abs(
vy_c(kc,j-1,i))) &
1752 vx_c_help = 0.5_dp*(
vx_c(kc,j,i)+
vx_c(kc,j,i-1))
1753 vy_c_help = 0.5_dp*(
vy_c(kc,j,i)+
vy_c(kc,j-1,i))
1755 lgs_b(kc) =
age_c(kc,j,i) + dtime_temp &
1757 ( min(vx_c_help, 0.0_dp) &
1760 +max(vx_c_help, 0.0_dp) &
1764 ( min(vy_c_help, 0.0_dp) &
1767 +max(vy_c_help, 0.0_dp) &
1777 lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
1779 lgs_a1(kc) = 1.0_dp+ci1(kc)*(ci2(kc)+ci2(kc-1))
1780 lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
1785 lgs_a0(kc) = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1))
1786 lgs_a1(kc) = 1.0_dp &
1787 +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
1788 -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
1789 lgs_a2(kc) = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
1793 adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
1795 lgs_a0(kc) = -max(adv_vert_help, 0.0_dp)
1796 lgs_a1(kc) = 1.0_dp &
1797 +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp)
1798 lgs_a2(kc) = min(adv_vert_help, 0.0_dp)
1804 lgs_b(kc) =
age_c(kc,j,i) + dtime_temp &
1806 ( (
vx_c(kc,j,i)-abs(
vx_c(kc,j,i))) &
1809 +(
vx_c(kc,j,i-1)+abs(
vx_c(kc,j,i-1))) &
1813 ( (
vy_c(kc,j,i)-abs(
vy_c(kc,j,i))) &
1816 +(
vy_c(kc,j-1,i)+abs(
vy_c(kc,j-1,i))) &
1822 vx_c_help = 0.5_dp*(
vx_c(kc,j,i)+
vx_c(kc,j,i-1))
1823 vy_c_help = 0.5_dp*(
vy_c(kc,j,i)+
vy_c(kc,j-1,i))
1825 lgs_b(kc) =
age_c(kc,j,i) + dtime_temp &
1827 ( min(vx_c_help, 0.0_dp) &
1830 +max(vx_c_help, 0.0_dp) &
1834 ( min(vy_c_help, 0.0_dp) &
1837 +max(vy_c_help, 0.0_dp) &
1846 if (
as_perp(j,i) >= zero)
then 1851 lgs_a0(kc) = -max(adv_vert_sg(kc-1), 0.0_dp)
1852 lgs_a1(kc) = 1.0_dp + max(adv_vert_sg(kc-1), 0.0_dp)
1857 lgs_b(kc) =
age_c(kc,j,i) + dtime_temp &
1859 ( (
vx_c(kc,j,i)-abs(
vx_c(kc,j,i))) &
1862 +(
vx_c(kc,j,i-1)+abs(
vx_c(kc,j,i-1))) &
1866 ( (
vy_c(kc,j,i)-abs(
vy_c(kc,j,i))) &
1869 +(
vy_c(kc,j-1,i)+abs(
vy_c(kc,j-1,i))) &
1875 vx_c_help = 0.5_dp*(
vx_c(kc,j,i)+
vx_c(kc,j,i-1))
1876 vy_c_help = 0.5_dp*(
vy_c(kc,j,i)+
vy_c(kc,j-1,i))
1878 lgs_b(kc) =
age_c(kc,j,i) + dtime_temp &
1880 ( min(vx_c_help, 0.0_dp) &
1883 +max(vx_c_help, 0.0_dp) &
1887 ( min(vy_c_help, 0.0_dp) &
1890 +max(vy_c_help, 0.0_dp) &
1900 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
1909 if (
age_c_neu(kc,j,i) < (age_min*year_sec)) &
1911 if (
age_c_neu(kc,j,i) > (age_max*year_sec)) &
1922 end subroutine calc_temp1
1928 subroutine calc_temp2(at1, at2_1, at2_2, at3_1, at3_2, &
1929 at4_1, at4_2, at5, at6, at7, atr1, alb1, &
1931 dtime_temp, dtt_2dxi, dtt_2deta, i, j)
1935 #if !defined(ALLOW_OPENAD) /* Normal */ 1939 #endif /* Normal vs. OpenAD */ 1943 #if !defined(ALLOW_OPENAD) /* Normal */ 1944 integer(i4b),
intent(in) :: i, j
1946 integer(i4b),
intent(inout) :: i, j
1947 #endif /* Normal vs. OpenAD */ 1948 real(dp),
intent(in) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
1949 at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
1950 at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
1951 ai1(0:kcmax), ai2(0:kcmax), &
1953 real(dp),
intent(in) :: dtime_temp, dtt_2dxi, dtt_2deta
1955 integer(i4b) :: kc, kt, kr
1956 real(dp) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), ct4(0:kcmax), &
1957 ct5(0:kcmax), ct6(0:kcmax), ct7(0:kcmax), ctr1, clb1
1958 real(dp) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), ct3_sg(0:kcmax), &
1959 ct4_sg(0:kcmax), adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
1960 real(dp) :: ci1(0:kcmax), ci2(0:kcmax)
1961 real(dp) :: temp_c_help(0:kcmax)
1962 real(dp) :: vx_c_help, vy_c_help
1963 real(dp) :: adv_vert_help
1964 real(dp) :: dtt_dxi, dtt_deta
1965 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
1966 lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
1967 lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
1968 lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
1969 lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
1970 real(dp),
parameter :: zero=0.0_dp
1974 if ((i == 0).or.(i == imax).or.(j == 0).or.(j == jmax))
then 1975 errormsg =
' >>> calc_temp2: Boundary points not allowed!' 1982 clb1 = alb1*
q_geo(j,i)
1987 ct1(kc) = at1(kc)/
h_c(j,i)*0.5_dp*(
vz_c(kc,j,i)+
vz_c(kc-1,j,i))
1991 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/
h_c(j,i)*
vz_c(kc,j,i)
1994 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/
h_c(j,i)*
vz_c(kc,j,i)
1997 #elif (ADV_VERT==2 || ADV_VERT==3) 2000 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/
h_c(j,i)*
vz_c(kc,j,i)
2007 ct2(kc) = ( at2_1(kc)*
dzm_dtau(j,i) &
2030 ct7(kc) = 2.0_dp*at7 &
2034 #if !defined(ALLOW_OPENAD) 2035 enh_c(kc,j,i), 0_i1b) &
2037 enh_c(kc,j,i), 0_i4b) &
2043 ci1(kc) = ai1(kc)/
h_c(j,i)
2050 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
2051 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
2052 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
2053 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
2054 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
2056 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
2057 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
2058 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
2059 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
2060 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
2062 #elif (ADV_VERT==2 || ADV_VERT==3) 2065 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
2066 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
2067 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
2068 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
2069 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
2075 temp_c_help(kc) = 0.5_dp*(
temp_c(kc,j,i)+
temp_c(kc+1,j,i))
2079 ci2(kc) = ai2(kc)/
h_c(j,i)
2083 dtt_dxi = 2.0_dp*dtt_2dxi
2084 dtt_deta = 2.0_dp*dtt_2deta
2091 lgs_a2(kr) = -1.0_dp
2099 lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
2101 lgs_b(kr) =
temp_r(kr,j,i)
2110 lgs_a2(kr) = -1.0_dp
2111 lgs_b(kr) = 2.0_dp*clb1
2123 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, krmax)
2142 lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2144 lgs_a1(kc) = 1.0_dp+ct5(kc)*(ct6(kc)+ct6(kc-1))
2145 lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2151 = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2155 +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2156 -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
2157 +ct5(kc)*(ct6(kc)+ct6(kc-1))
2159 = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
2164 adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
2167 = -max(adv_vert_help, 0.0_dp) &
2171 +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp) &
2172 +ct5(kc)*(ct6(kc)+ct6(kc-1))
2174 = min(adv_vert_help, 0.0_dp) &
2181 lgs_b(kc) =
temp_c(kc,j,i) + ct7(kc) &
2183 ( (
vx_c(kc,j,i)-abs(
vx_c(kc,j,i))) &
2186 +(
vx_c(kc,j,i-1)+abs(
vx_c(kc,j,i-1))) &
2190 ( (
vy_c(kc,j,i)-abs(
vy_c(kc,j,i))) &
2193 +(
vy_c(kc,j-1,i)+abs(
vy_c(kc,j-1,i))) &
2199 vx_c_help = 0.5_dp*(
vx_c(kc,j,i)+
vx_c(kc,j,i-1))
2200 vy_c_help = 0.5_dp*(
vy_c(kc,j,i)+
vy_c(kc,j-1,i))
2202 lgs_b(kc) =
temp_c(kc,j,i) + ct7(kc) &
2204 ( min(vx_c_help, 0.0_dp) &
2207 +max(vx_c_help, 0.0_dp) &
2211 ( min(vy_c_help, 0.0_dp) &
2214 +max(vy_c_help, 0.0_dp) &
2229 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
2251 lgs_a1(kc) = 1.0_dp - min(adv_vert_sg(kc), 0.0_dp)
2252 lgs_a2(kc) = min(adv_vert_sg(kc), 0.0_dp)
2256 lgs_b(kc) =
age_c(kc,j,i) + dtime_temp &
2258 ( (
vx_c(kc,j,i)-abs(
vx_c(kc,j,i))) &
2261 +(
vx_c(kc,j,i-1)+abs(
vx_c(kc,j,i-1))) &
2265 ( (
vy_c(kc,j,i)-abs(
vy_c(kc,j,i))) &
2268 +(
vy_c(kc,j-1,i)+abs(
vy_c(kc,j-1,i))) &
2274 vx_c_help = 0.5_dp*(
vx_c(kc,j,i)+
vx_c(kc,j,i-1))
2275 vy_c_help = 0.5_dp*(
vy_c(kc,j,i)+
vy_c(kc,j-1,i))
2277 lgs_b(kc) =
age_c(kc,j,i) + dtime_temp &
2279 ( min(vx_c_help, 0.0_dp) &
2282 +max(vx_c_help, 0.0_dp) &
2286 ( min(vy_c_help, 0.0_dp) &
2289 +max(vy_c_help, 0.0_dp) &
2299 lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2301 lgs_a1(kc) = 1.0_dp+ci1(kc)*(ci2(kc)+ci2(kc-1))
2302 lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2307 lgs_a0(kc) = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1))
2308 lgs_a1(kc) = 1.0_dp &
2309 +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2310 -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
2311 lgs_a2(kc) = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
2315 adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
2317 lgs_a0(kc) = -max(adv_vert_help, 0.0_dp)
2318 lgs_a1(kc) = 1.0_dp &
2319 +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp)
2320 lgs_a2(kc) = min(adv_vert_help, 0.0_dp)
2326 lgs_b(kc) =
age_c(kc,j,i) + dtime_temp &
2328 ( (
vx_c(kc,j,i)-abs(
vx_c(kc,j,i))) &
2331 +(
vx_c(kc,j,i-1)+abs(
vx_c(kc,j,i-1))) &
2335 ( (
vy_c(kc,j,i)-abs(
vy_c(kc,j,i))) &
2338 +(
vy_c(kc,j-1,i)+abs(
vy_c(kc,j-1,i))) &
2344 vx_c_help = 0.5_dp*(
vx_c(kc,j,i)+
vx_c(kc,j,i-1))
2345 vy_c_help = 0.5_dp*(
vy_c(kc,j,i)+
vy_c(kc,j-1,i))
2347 lgs_b(kc) =
age_c(kc,j,i) + dtime_temp &
2349 ( min(vx_c_help, 0.0_dp) &
2352 +max(vx_c_help, 0.0_dp) &
2356 ( min(vy_c_help, 0.0_dp) &
2359 +max(vy_c_help, 0.0_dp) &
2368 if (
as_perp(j,i) >= zero)
then 2373 lgs_a0(kc) = -max(adv_vert_sg(kc-1), 0.0_dp)
2374 lgs_a1(kc) = 1.0_dp + max(adv_vert_sg(kc-1), 0.0_dp)
2379 lgs_b(kc) =
age_c(kc,j,i) + dtime_temp &
2381 ( (
vx_c(kc,j,i)-abs(
vx_c(kc,j,i))) &
2384 +(
vx_c(kc,j,i-1)+abs(
vx_c(kc,j,i-1))) &
2388 ( (
vy_c(kc,j,i)-abs(
vy_c(kc,j,i))) &
2391 +(
vy_c(kc,j-1,i)+abs(
vy_c(kc,j-1,i))) &
2397 vx_c_help = 0.5_dp*(
vx_c(kc,j,i)+
vx_c(kc,j,i-1))
2398 vy_c_help = 0.5_dp*(
vy_c(kc,j,i)+
vy_c(kc,j-1,i))
2400 lgs_b(kc) =
age_c(kc,j,i) + dtime_temp &
2402 ( min(vx_c_help, 0.0_dp) &
2405 +max(vx_c_help, 0.0_dp) &
2409 ( min(vy_c_help, 0.0_dp) &
2412 +max(vy_c_help, 0.0_dp) &
2422 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
2431 if (
age_c_neu(kc,j,i) < (age_min*year_sec)) &
2433 if (
age_c_neu(kc,j,i) > (age_max*year_sec)) &
2450 subroutine calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
2451 at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, &
2452 aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
2453 ai1, ai2, ai3, dzeta_t, &
2454 dtime_temp, dtt_2dxi, dtt_2deta, i, j)
2458 #if !defined(ALLOW_OPENAD) /* Normal */ 2462 #endif /* Normal vs. OpenAD */ 2466 #if !defined(ALLOW_OPENAD) /* Normal */ 2467 integer(i4b),
intent(in) :: i, j
2469 integer(i4b),
intent(inout) :: i, j
2470 #endif /* Normal vs. OpenAD */ 2471 real(dp),
intent(in) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
2472 at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
2473 at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
2474 ai1(0:kcmax), ai2(0:kcmax), ai3, &
2475 atr1, am1, am2, alb1
2476 real(dp),
intent(in) :: aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld
2477 real(dp),
intent(in) :: dzeta_t, dtime_temp, dtt_2dxi, dtt_2deta
2479 integer(i4b) :: kc, kt, kr
2480 real(dp) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), ct4(0:kcmax), &
2481 ct5(0:kcmax), ct6(0:kcmax), ct7(0:kcmax), ctr1, cm1, cm2, &
2483 real(dp) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), ct3_sg(0:kcmax), &
2484 ct4_sg(0:kcmax), adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
2485 real(dp) :: ci1(0:kcmax), ci2(0:kcmax), ci3
2486 real(dp) :: cw1(0:ktmax), cw2(0:ktmax), cw3(0:ktmax), cw4(0:ktmax), &
2487 cw5, cw7(0:ktmax), cw8, cw9(0:ktmax)
2488 real(dp) :: cw1_sg(0:ktmax), cw2_sg(0:ktmax), cw3_sg(0:ktmax), &
2489 cw4_sg(0:ktmax), adv_vert_w_sg(0:ktmax), abs_adv_vert_w_sg(0:ktmax)
2490 real(dp) :: sigma_c_help(0:kcmax), sigma_t_help(0:ktmax), &
2491 temp_c_help(0:kcmax)
2492 real(dp) :: vx_c_help, vy_c_help, vx_t_help, vy_t_help
2493 real(dp) :: adv_vert_help, adv_vert_w_help
2494 real(dp) :: dtt_dxi, dtt_deta
2495 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
2496 lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
2497 lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
2498 lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
2499 lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
2500 real(dp),
parameter :: zero=0.0_dp
2504 if ((i == 0).or.(i == imax).or.(j == 0).or.(j == jmax))
then 2505 errormsg =
' >>> calc_temp3: Boundary points not allowed!' 2513 clb1 = alb1*
q_geo(j,i)
2522 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/
h_c_neu(j,i)*
vz_c(kc,j,i)
2525 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/
h_c_neu(j,i)*
vz_c(kc,j,i)
2528 #elif (ADV_VERT==2 || ADV_VERT==3) 2531 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/
h_c_neu(j,i)*
vz_c(kc,j,i)
2538 ct2(kc) = ( at2_1(kc)*
dzm_dtau(j,i) &
2561 *
creep(sigma_c_help(kc)) &
2562 *sigma_c_help(kc)*sigma_c_help(kc)
2565 ct7(kc) = 2.0_dp*at7 &
2569 #if !defined(ALLOW_OPENAD) 2570 enh_c(kc,j,i), 0_i1b) &
2572 enh_c(kc,j,i), 0_i4b) &
2578 ci1(kc) = ai1(kc)/
h_c_neu(j,i)
2585 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
2586 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
2587 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
2588 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
2589 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
2591 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
2592 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
2593 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
2594 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
2595 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
2597 #elif (ADV_VERT==2 || ADV_VERT==3) 2600 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
2601 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
2602 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
2603 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
2604 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
2610 temp_c_help(kc) = 0.5_dp*(
temp_c(kc,j,i)+
temp_c(kc+1,j,i))
2614 ci2(kc) = ai2(kc)/
h_c_neu(j,i)
2637 #elif (ADV_VERT==2 || ADV_VERT==3) 2651 -0.5_dp*(
vz_t(kt,j,i)+
vz_t(kt-1,j,i)) )
2675 *
creep(sigma_t_help(kt)) &
2676 *sigma_t_help(kt)*sigma_t_help(kt)
2679 cw7(kt) = 2.0_dp*aw7 &
2683 #if !defined(ALLOW_OPENAD) 2684 enh_t(kt,j,i), 0_i1b) &
2686 enh_t(kt,j,i), 0_i4b) &
2697 cw2_sg(kt) = 0.5_dp*(cw2(kt)+cw2(kt+1))
2698 cw3_sg(kt) = 0.5_dp*(cw3(kt)+cw3(kt+1))
2699 cw4_sg(kt) = 0.5_dp*(cw4(kt)+cw4(kt+1))
2700 adv_vert_w_sg(kt) = cw1_sg(kt)-cw2_sg(kt)-cw3_sg(kt)-cw4_sg(kt)
2701 abs_adv_vert_w_sg(kt) = abs(adv_vert_w_sg(kt))
2703 cw2_sg(kt) = 0.5_dp*(cw2(kt)+cw2(kt+1))
2704 cw3_sg(kt) = 0.5_dp*(cw3(kt)+cw3(kt+1))
2705 cw4_sg(kt) = 0.5_dp*(cw4(kt)+cw4(kt+1))
2706 adv_vert_w_sg(kt) = cw1_sg(kt)-cw2_sg(kt)-cw3_sg(kt)-cw4_sg(kt)
2707 abs_adv_vert_w_sg(kt) = abs(adv_vert_w_sg(kt))
2709 #elif (ADV_VERT==2 || ADV_VERT==3) 2712 cw2_sg(kt) = 0.5_dp*(cw2(kt)+cw2(kt+1))
2713 cw3_sg(kt) = 0.5_dp*(cw3(kt)+cw3(kt+1))
2714 cw4_sg(kt) = 0.5_dp*(cw4(kt)+cw4(kt+1))
2715 adv_vert_w_sg(kt) = cw1_sg(kt)-cw2_sg(kt)-cw3_sg(kt)-cw4_sg(kt)
2716 abs_adv_vert_w_sg(kt) = abs(adv_vert_w_sg(kt))
2722 dtt_dxi = 2.0_dp*dtt_2dxi
2723 dtt_deta = 2.0_dp*dtt_2deta
2730 lgs_a2(kr) = -1.0_dp
2738 lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
2740 lgs_b(kr) =
temp_r(kr,j,i)
2749 lgs_a2(kr) = -1.0_dp
2750 lgs_b(kr) = 2.0_dp*clb1
2762 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, krmax)
2775 lgs_a2(kt) = -1.0_dp
2782 lgs_a0(kt) = -0.5_dp*(cw1(kt)-cw2(kt)-cw3(kt)-cw4(kt)) - cw5
2783 lgs_a1(kt) = 1.0_dp + 2.0_dp*cw5
2784 lgs_a2(kt) = 0.5_dp*(cw1(kt)-cw2(kt)-cw3(kt)-cw4(kt)) - cw5
2789 = -0.5_dp*(adv_vert_w_sg(kt-1)+abs_adv_vert_w_sg(kt-1)) &
2793 +0.5_dp*(adv_vert_w_sg(kt-1)+abs_adv_vert_w_sg(kt-1)) &
2794 -0.5_dp*(adv_vert_w_sg(kt) -abs_adv_vert_w_sg(kt) ) &
2797 = 0.5_dp*(adv_vert_w_sg(kt) -abs_adv_vert_w_sg(kt) ) &
2802 adv_vert_w_help = 0.5_dp*(adv_vert_w_sg(kt)+adv_vert_w_sg(kt-1))
2805 = -max(adv_vert_w_help, 0.0_dp) &
2809 +max(adv_vert_w_help, 0.0_dp)-min(adv_vert_w_help, 0.0_dp) &
2812 = min(adv_vert_w_help, 0.0_dp) &
2819 lgs_b(kt) =
omega_t(kt,j,i) + cw7(kt) + cw8 + cw9(kt) &
2821 ( (
vx_t(kt,j,i)-abs(
vx_t(kt,j,i))) &
2824 +(
vx_t(kt,j,i-1)+abs(
vx_t(kt,j,i-1))) &
2828 ( (
vy_t(kt,j,i)-abs(
vy_t(kt,j,i))) &
2831 +(
vy_t(kt,j-1,i)+abs(
vy_t(kt,j-1,i))) &
2837 vx_t_help = 0.5_dp*(
vx_t(kt,j,i)+
vx_t(kt,j,i-1))
2838 vy_t_help = 0.5_dp*(
vy_t(kt,j,i)+
vy_t(kt,j-1,i))
2840 lgs_b(kt) =
omega_t(kt,j,i) + cw7(kt) + cw8 + cw9(kt) &
2842 ( min(vx_t_help, 0.0_dp) &
2845 +max(vx_t_help, 0.0_dp) &
2849 ( min(vy_t_help, 0.0_dp) &
2852 +max(vy_t_help, 0.0_dp) &
2862 #if (!defined(CTS_MELTING_FREEZING) || CTS_MELTING_FREEZING==1) 2864 if (
am_perp(j,i) >= zero)
then 2869 lgs_a0(kt) = -1.0_dp
2874 #elif (CTS_MELTING_FREEZING==2) 2882 errormsg =
' >>> calc_temp3: CTS_MELTING_FREEZING must be either 1 or 2!' 2889 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, ktmax)
2897 if (lgs_x(kt) < zero)
then 2914 #if (!defined(CTS_MELTING_FREEZING) || CTS_MELTING_FREEZING==1) 2916 if (
am_perp(j,i) >= zero)
then 2923 #elif (CTS_MELTING_FREEZING==2) 2929 errormsg =
' >>> calc_temp3: CTS_MELTING_FREEZING must be either 1 or 2!' 2936 lgs_a2(kc) = -1.0_dp
2937 lgs_b(kc) = -cm1-cm2
2943 lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2945 lgs_a1(kc) = 1.0_dp+ct5(kc)*(ct6(kc)+ct6(kc-1))
2946 lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
2952 = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2956 +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
2957 -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
2958 +ct5(kc)*(ct6(kc)+ct6(kc-1))
2960 = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
2965 adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
2968 = -max(adv_vert_help, 0.0_dp) &
2972 +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp) &
2973 +ct5(kc)*(ct6(kc)+ct6(kc-1))
2975 = min(adv_vert_help, 0.0_dp) &
2982 lgs_b(kc) =
temp_c(kc,j,i) + ct7(kc) &
2984 ( (
vx_c(kc,j,i)-abs(
vx_c(kc,j,i))) &
2987 +(
vx_c(kc,j,i-1)+abs(
vx_c(kc,j,i-1))) &
2991 ( (
vy_c(kc,j,i)-abs(
vy_c(kc,j,i))) &
2994 +(
vy_c(kc,j-1,i)+abs(
vy_c(kc,j-1,i))) &
3000 vx_c_help = 0.5_dp*(
vx_c(kc,j,i)+
vx_c(kc,j,i-1))
3001 vy_c_help = 0.5_dp*(
vy_c(kc,j,i)+
vy_c(kc,j-1,i))
3003 lgs_b(kc) =
temp_c(kc,j,i) + ct7(kc) &
3005 ( min(vx_c_help, 0.0_dp) &
3008 +max(vx_c_help, 0.0_dp) &
3012 ( min(vy_c_help, 0.0_dp) &
3015 +max(vy_c_help, 0.0_dp) &
3030 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
3042 lgs_a1(kt) = 1.0_dp - min(adv_vert_w_sg(kt), 0.0_dp)
3043 lgs_a2(kt) = min(adv_vert_w_sg(kt), 0.0_dp)
3047 lgs_b(kt) =
age_t(kt,j,i) + dtime_temp &
3049 ( (
vx_t(kt,j,i)-abs(
vx_t(kt,j,i))) &
3052 +(
vx_t(kt,j,i-1)+abs(
vx_t(kt,j,i-1))) &
3056 ( (
vy_t(kt,j,i)-abs(
vy_t(kt,j,i))) &
3059 +(
vy_t(kt,j-1,i)+abs(
vy_t(kt,j-1,i))) &
3065 vx_t_help = 0.5_dp*(
vx_t(kt,j,i)+
vx_t(kt,j,i-1))
3066 vy_t_help = 0.5_dp*(
vy_t(kt,j,i)+
vy_t(kt,j-1,i))
3068 lgs_b(kt) =
age_t(kt,j,i) + dtime_temp &
3070 ( min(vx_t_help, 0.0_dp) &
3073 +max(vx_t_help, 0.0_dp) &
3077 ( min(vy_t_help, 0.0_dp) &
3080 +max(vy_t_help, 0.0_dp) &
3090 lgs_a0(kt) = -0.5_dp*(cw1(kt)-cw2(kt)-cw3(kt)-cw4(kt)) - ci3
3091 lgs_a1(kt) = 1.0_dp + 2.0_dp*ci3
3092 lgs_a2(kt) = 0.5_dp*(cw1(kt)-cw2(kt)-cw3(kt)-cw4(kt)) - ci3
3096 lgs_a0(kt) = -0.5_dp*(adv_vert_w_sg(kt-1)+abs_adv_vert_w_sg(kt-1))
3097 lgs_a1(kt) = 1.0_dp &
3098 +0.5_dp*(adv_vert_w_sg(kt-1)+abs_adv_vert_w_sg(kt-1)) &
3099 -0.5_dp*(adv_vert_w_sg(kt) -abs_adv_vert_w_sg(kt) )
3100 lgs_a2(kt) = 0.5_dp*(adv_vert_w_sg(kt) -abs_adv_vert_w_sg(kt) )
3104 adv_vert_w_help = 0.5_dp*(adv_vert_w_sg(kt)+adv_vert_w_sg(kt-1))
3106 lgs_a0(kt) = -max(adv_vert_w_help, 0.0_dp)
3107 lgs_a1(kt) = 1.0_dp &
3108 +max(adv_vert_w_help, 0.0_dp)-min(adv_vert_w_help, 0.0_dp)
3109 lgs_a2(kt) = min(adv_vert_w_help, 0.0_dp)
3115 lgs_b(kt) =
age_t(kt,j,i) + dtime_temp &
3117 ( (
vx_t(kt,j,i)-abs(
vx_t(kt,j,i))) &
3120 +(
vx_t(kt,j,i-1)+abs(
vx_t(kt,j,i-1))) &
3124 ( (
vy_t(kt,j,i)-abs(
vy_t(kt,j,i))) &
3127 +(
vy_t(kt,j-1,i)+abs(
vy_t(kt,j-1,i))) &
3133 vx_t_help = 0.5_dp*(
vx_t(kt,j,i)+
vx_t(kt,j,i-1))
3134 vy_t_help = 0.5_dp*(
vy_t(kt,j,i)+
vy_t(kt,j-1,i))
3136 lgs_b(kt) =
age_t(kt,j,i) + dtime_temp &
3138 ( min(vx_t_help, 0.0_dp) &
3141 +max(vx_t_help, 0.0_dp) &
3145 ( min(vy_t_help, 0.0_dp) &
3148 +max(vy_t_help, 0.0_dp) &
3161 lgs_a0(kt) = -0.5_dp*(cw1(kt)-cw2(kt)-cw3(kt)-cw4(kt)) - ci3
3162 lgs_a1(kt) = 1.0_dp + 2.0_dp*ci3
3163 lgs_a2(kt) = 0.5_dp*(cw1(kt)-cw2(kt)-cw3(kt)-cw4(kt)) - ci3
3167 lgs_b(kt) =
age_t(kt,j,i) + dtime_temp &
3169 ( (
vx_t(kt,j,i)-abs(
vx_t(kt,j,i))) &
3172 +(
vx_t(kt,j,i-1)+abs(
vx_t(kt,j,i-1))) &
3176 ( (
vy_t(kt,j,i)-abs(
vy_t(kt,j,i))) &
3179 +(
vy_t(kt,j-1,i)+abs(
vy_t(kt,j-1,i))) &
3185 vx_t_help = 0.5_dp*(
vx_t(kt,j,i)+
vx_t(kt,j,i-1))
3186 vy_t_help = 0.5_dp*(
vy_t(kt,j,i)+
vy_t(kt,j-1,i))
3188 lgs_b(kt) =
age_t(kt,j,i) + dtime_temp &
3190 ( min(vx_t_help, 0.0_dp) &
3193 +max(vx_t_help, 0.0_dp) &
3197 ( min(vy_t_help, 0.0_dp) &
3200 +max(vy_t_help, 0.0_dp) &
3206 #elif (ADV_VERT==2 || ADV_VERT==3) 3211 if (adv_vert_sg(kc) <= zero)
then 3213 lgs_a0(ktmax+kc) = 0.0_dp
3214 lgs_a1(ktmax+kc) = 1.0_dp - adv_vert_sg(kc)
3215 lgs_a2(ktmax+kc) = adv_vert_sg(kc)
3219 lgs_b(ktmax+kc) =
age_c(kc,j,i) + dtime_temp &
3221 ( (
vx_c(kc,j,i)-abs(
vx_c(kc,j,i))) &
3224 +(
vx_c(kc,j,i-1)+abs(
vx_c(kc,j,i-1))) &
3228 ( (
vy_c(kc,j,i)-abs(
vy_c(kc,j,i))) &
3231 +(
vy_c(kc,j-1,i)+abs(
vy_c(kc,j-1,i))) &
3237 vx_c_help = 0.5_dp*(
vx_c(kc,j,i)+
vx_c(kc,j,i-1))
3238 vy_c_help = 0.5_dp*(
vy_c(kc,j,i)+
vy_c(kc,j-1,i))
3240 lgs_b(ktmax+kc) =
age_c(kc,j,i) + dtime_temp &
3242 ( min(vx_c_help, 0.0_dp) &
3245 +max(vx_c_help, 0.0_dp) &
3249 ( min(vy_c_help, 0.0_dp) &
3252 +max(vy_c_help, 0.0_dp) &
3258 else if (adv_vert_w_sg(kt-1) >= zero)
then 3260 lgs_a0(kt) = -adv_vert_w_sg(kt-1)
3261 lgs_a1(kt) = 1.0_dp + adv_vert_w_sg(kt-1)
3266 lgs_b(kt) =
age_t(kt,j,i) + dtime_temp &
3268 ( (
vx_t(kt,j,i)-abs(
vx_t(kt,j,i))) &
3271 +(
vx_t(kt,j,i-1)+abs(
vx_t(kt,j,i-1))) &
3275 ( (
vy_t(kt,j,i)-abs(
vy_t(kt,j,i))) &
3278 +(
vy_t(kt,j-1,i)+abs(
vy_t(kt,j-1,i))) &
3284 vx_t_help = 0.5_dp*(
vx_t(kt,j,i)+
vx_t(kt,j,i-1))
3285 vy_t_help = 0.5_dp*(
vy_t(kt,j,i)+
vy_t(kt,j-1,i))
3287 lgs_b(kt) =
age_t(kt,j,i) + dtime_temp &
3289 ( min(vx_t_help, 0.0_dp) &
3292 +max(vx_t_help, 0.0_dp) &
3296 ( min(vy_t_help, 0.0_dp) &
3299 +max(vy_t_help, 0.0_dp) &
3307 lgs_a0(kt) = -0.5_dp
3309 lgs_a2(kt) = -0.5_dp
3321 lgs_a0(ktmax+kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
3323 lgs_a1(ktmax+kc) = 1.0_dp+ci1(kc)*(ci2(kc)+ci2(kc-1))
3324 lgs_a2(ktmax+kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
3329 lgs_a0(ktmax+kc) = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1))
3330 lgs_a1(ktmax+kc) = 1.0_dp &
3331 +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
3332 -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
3333 lgs_a2(ktmax+kc) = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
3337 adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
3339 lgs_a0(ktmax+kc) = -max(adv_vert_help, 0.0_dp)
3340 lgs_a1(ktmax+kc) = 1.0_dp &
3341 +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp)
3342 lgs_a2(ktmax+kc) = min(adv_vert_help, 0.0_dp)
3348 lgs_b(ktmax+kc) =
age_c(kc,j,i) + dtime_temp &
3350 ( (
vx_c(kc,j,i)-abs(
vx_c(kc,j,i))) &
3353 +(
vx_c(kc,j,i-1)+abs(
vx_c(kc,j,i-1))) &
3357 ( (
vy_c(kc,j,i)-abs(
vy_c(kc,j,i))) &
3360 +(
vy_c(kc,j-1,i)+abs(
vy_c(kc,j-1,i))) &
3366 vx_c_help = 0.5_dp*(
vx_c(kc,j,i)+
vx_c(kc,j,i-1))
3367 vy_c_help = 0.5_dp*(
vy_c(kc,j,i)+
vy_c(kc,j-1,i))
3369 lgs_b(ktmax+kc) =
age_c(kc,j,i) + dtime_temp &
3371 ( min(vx_c_help, 0.0_dp) &
3374 +max(vx_c_help, 0.0_dp) &
3378 ( min(vy_c_help, 0.0_dp) &
3381 +max(vy_c_help, 0.0_dp) &
3390 if (
as_perp(j,i) >= zero)
then 3391 lgs_a0(ktmax+kc) = 0.0_dp
3392 lgs_a1(ktmax+kc) = 1.0_dp
3393 lgs_b(ktmax+kc) = 0.0_dp
3395 lgs_a0(ktmax+kc) = -max(adv_vert_sg(kc-1), 0.0_dp)
3396 lgs_a1(ktmax+kc) = 1.0_dp + max(adv_vert_sg(kc-1), 0.0_dp)
3401 lgs_b(ktmax+kc) =
age_c(kc,j,i) + dtime_temp &
3403 ( (
vx_c(kc,j,i)-abs(
vx_c(kc,j,i))) &
3406 +(
vx_c(kc,j,i-1)+abs(
vx_c(kc,j,i-1))) &
3410 ( (
vy_c(kc,j,i)-abs(
vy_c(kc,j,i))) &
3413 +(
vy_c(kc,j-1,i)+abs(
vy_c(kc,j-1,i))) &
3419 vx_c_help = 0.5_dp*(
vx_c(kc,j,i)+
vx_c(kc,j,i-1))
3420 vy_c_help = 0.5_dp*(
vy_c(kc,j,i)+
vy_c(kc,j-1,i))
3422 lgs_b(ktmax+kc) =
age_c(kc,j,i) + dtime_temp &
3424 ( min(vx_c_help, 0.0_dp) &
3427 +max(vx_c_help, 0.0_dp) &
3431 ( min(vy_c_help, 0.0_dp) &
3434 +max(vy_c_help, 0.0_dp) &
3444 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax+ktmax)
3453 if (
age_t_neu(kt,j,i) < (age_min*year_sec)) &
3455 if (
age_t_neu(kt,j,i) > (age_max*year_sec)) &
3464 if (
age_c_neu(kc,j,i) < (age_min*year_sec)) &
3466 if (
age_c_neu(kc,j,i) > (age_max*year_sec)) &
3478 #if !defined(ALLOW_OPENAD) /* Normal */ 3482 #endif /* Normal vs. OpenAD */ 3486 #if !defined(ALLOW_OPENAD) /* Normal */ 3487 integer(i4b),
intent(in) :: i, j
3489 integer(i4b),
intent(inout) :: i, j
3490 #endif /* Normal vs. OpenAD */ 3491 real(dp),
intent(in) :: atr1, alb1
3493 integer(i4b) :: kc, kt, kr
3494 real(dp) :: ctr1, clb1
3495 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
3496 lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
3497 lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
3498 lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
3499 lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
3504 clb1 = alb1*
q_geo(j,i)
3509 #if defined(ALLOW_OPENAD) /* OpenAD */ 3513 lgs_a2(kr) = -1.0_dp
3521 lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
3523 lgs_b(kr) =
temp_r(kr,j,i)
3532 lgs_a2(kr) = -1.0_dp
3533 lgs_b(kr) = 2.0_dp*clb1
3545 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, krmax)
3575 at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, &
3576 aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
3577 ai1, ai2, ai3, dzeta_t, &
3578 dtime_temp, dtt_2dxi, dtt_2deta, dtime_temp_inv, &
3583 #if !defined(ALLOW_OPENAD) /* Normal */ 3584 integer(i4b),
intent(in) :: i, j
3586 integer(i4b),
intent(inout) :: i, j
3587 #endif /* Normal vs. OpenAD */ 3588 real(dp),
intent(in) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
3589 at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
3590 at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
3591 ai1(0:kcmax), ai2(0:kcmax), ai3, &
3592 atr1, am1, am2, alb1
3593 real(dp),
intent(in) :: aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld
3594 real(dp),
intent(in) :: dzeta_t
3595 real(dp),
intent(in) :: dtime_temp, dtime_temp_inv, dtt_2dxi, dtt_2deta
3597 real(dp) :: zm_shift
3598 real(dp) :: difftemp_a, difftemp_b, interpol
3606 if (difftemp_a <= 0.0_dp)
return 3610 do while (difftemp_a > 0.0_dp)
3626 call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
3627 at4_1, at4_2, at5, at6, at7, atr1, &
3629 aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
3630 ai1, ai2, ai3, dzeta_t, &
3631 dtime_temp, dtt_2dxi, dtt_2deta, i, j)
3633 difftemp_b = difftemp_a
3642 interpol = difftemp_a/(difftemp_b-difftemp_a)*zm_shift
3654 call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
3655 at4_1, at4_2, at5, at6, at7, atr1, &
3657 aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
3658 ai1, ai2, ai3, dzeta_t, &
3659 dtime_temp, dtt_2dxi, dtt_2deta, i, j)
3666 subroutine shift_cts_downward(at1, at2_1, at2_2, at3_1, at3_2, &
3667 at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, &
3668 aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
3669 ai1, ai2, ai3, dzeta_t, &
3670 dtime_temp, dtt_2dxi, dtt_2deta, dtime_temp_inv, &
3675 #if !defined(ALLOW_OPENAD) /* Normal */ 3676 integer(i4b),
intent(in) :: i, j
3678 integer(i4b),
intent(inout) :: i, j
3679 #endif /* Normal vs. OpenAD */ 3680 real(dp),
intent(in) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
3681 at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
3682 at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
3683 ai1(0:kcmax), ai2(0:kcmax), ai3, &
3684 atr1, am1, am2, alb1
3685 real(dp),
intent(in) :: aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld
3686 real(dp),
intent(in) :: dzeta_t
3687 real(dp),
intent(in) :: dtt_2dxi, dtt_2deta, dtime_temp, dtime_temp_inv
3689 real(dp) :: zm_shift
3690 real(dp) :: difftemp_a, difftemp_b, interpol
3698 if (difftemp_a >= 0.0_dp)
return 3702 do while (difftemp_a < 0.0_dp)
3710 zm_shift = (
zm_neu(j,i)+zm_shift)-(
zb(j,i)+0.001_dp)
3722 call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
3723 at4_1, at4_2, at5, at6, at7, atr1, &
3725 aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
3726 ai1, ai2, ai3, dzeta_t, &
3727 dtime_temp, dtt_2dxi, dtt_2deta, i, j)
3729 difftemp_b = difftemp_a
3732 if (difftemp_a >= 0.0_dp)
then 3738 interpol = difftemp_a/(difftemp_a-difftemp_b)*zm_shift
3750 call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
3751 at4_1, at4_2, at5, at6, at7, atr1, &
3753 aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
3754 ai1, ai2, ai3, dzeta_t, &
3755 dtime_temp, dtt_2dxi, dtt_2deta, i, j)
3770 call calc_temp2(at1, at2_1, at2_2, at3_1, at3_2, &
3771 at4_1, at4_2, at5, at6, at7, atr1, alb1, &
3773 dtime_temp, dtt_2dxi, dtt_2deta, i, j)
3791 call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
3792 at4_1, at4_2, at5, at6, at7, atr1, &
3794 aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
3795 ai1, ai2, ai3, dzeta_t, &
3796 dtime_temp, dtt_2dxi, dtt_2deta, i, j)
3798 difftemp_b = difftemp_a
3807 interpol = difftemp_a/(difftemp_a-difftemp_b)*zm_shift
3819 call calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, &
3820 at4_1, at4_2, at5, at6, at7, atr1, &
3822 aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, &
3823 ai1, ai2, ai3, dzeta_t, &
3824 dtime_temp, dtt_2dxi, dtt_2deta, i, j)
3826 end subroutine shift_cts_downward
3831 subroutine calc_temp_ssa(at1, at2_1, at2_2, at3_1, at3_2, &
3832 at4_1, at4_2, at5, at6, at7, atr1, alb1, &
3834 dtime_temp, dtt_2dxi, dtt_2deta, i, j)
3838 #if !defined(ALLOW_OPENAD) /* Normal */ 3842 #endif /* Normal vs. OpenAD */ 3846 #if !defined(ALLOW_OPENAD) /* Normal */ 3847 integer(i4b),
intent(in) :: i, j
3849 integer(i4b),
intent(inout) :: i, j
3850 #endif /* Normal vs. OpenAD */ 3851 real(dp),
intent(in) :: at1(0:kcmax), at2_1(0:kcmax), at2_2(0:kcmax), &
3852 at3_1(0:kcmax), at3_2(0:kcmax), at4_1(0:kcmax), &
3853 at4_2(0:kcmax), at5(0:kcmax), at6(0:kcmax), at7, &
3854 ai1(0:kcmax), ai2(0:kcmax), &
3856 real(dp),
intent(in) :: dtime_temp, dtt_2dxi, dtt_2deta
3858 integer(i4b) :: kc, kt, kr
3859 real(dp) :: ct1(0:kcmax), ct2(0:kcmax), ct3(0:kcmax), ct4(0:kcmax), &
3860 ct5(0:kcmax), ct6(0:kcmax), ct7(0:kcmax), ctr1, clb1
3861 real(dp) :: ct1_sg(0:kcmax), ct2_sg(0:kcmax), ct3_sg(0:kcmax), &
3862 ct4_sg(0:kcmax), adv_vert_sg(0:kcmax), abs_adv_vert_sg(0:kcmax)
3863 real(dp) :: ci1(0:kcmax), ci2(0:kcmax)
3864 real(dp) :: temp_c_help(0:kcmax)
3865 real(dp) :: vx_c_help, vy_c_help
3866 real(dp) :: adv_vert_help
3867 real(dp) :: dtt_dxi, dtt_deta
3868 real(dp) :: lgs_a0(0:kcmax+ktmax+krmax+imax+jmax), &
3869 lgs_a1(0:kcmax+ktmax+krmax+imax+jmax), &
3870 lgs_a2(0:kcmax+ktmax+krmax+imax+jmax), &
3871 lgs_x(0:kcmax+ktmax+krmax+imax+jmax), &
3872 lgs_b(0:kcmax+ktmax+krmax+imax+jmax)
3873 real(dp),
parameter :: zero=0.0_dp
3877 if ((i == 0).or.(i == imax).or.(j == 0).or.(j == jmax))
then 3878 errormsg =
' >>> calc_temp_ssa: Boundary points not allowed!' 3885 clb1 = alb1*
q_geo(j,i)
3890 ct1(kc) = at1(kc)/
h_c(j,i)*0.5_dp*(
vz_c(kc,j,i)+
vz_c(kc-1,j,i))
3894 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/
h_c(j,i)*
vz_c(kc,j,i)
3897 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/
h_c(j,i)*
vz_c(kc,j,i)
3900 #elif (ADV_VERT==2 || ADV_VERT==3) 3903 ct1_sg(kc) = 0.5_dp*(at1(kc)+at1(kc+1))/
h_c(j,i)*
vz_c(kc,j,i)
3910 ct2(kc) = ( at2_1(kc)*
dzm_dtau(j,i) &
3921 ct7(kc) = 2.0_dp*at7 &
3925 #if !defined(ALLOW_OPENAD) 3926 enh_c(kc,j,i), 0_i1b) &
3928 enh_c(kc,j,i), 0_i4b) &
3931 ci1(kc) = ai1(kc)/
h_c(j,i)
3938 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
3939 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
3940 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
3941 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
3942 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
3944 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
3945 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
3946 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
3947 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
3948 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
3950 #elif (ADV_VERT==2 || ADV_VERT==3) 3953 ct2_sg(kc) = 0.5_dp*(ct2(kc)+ct2(kc+1))
3954 ct3_sg(kc) = 0.5_dp*(ct3(kc)+ct3(kc+1))
3955 ct4_sg(kc) = 0.5_dp*(ct4(kc)+ct4(kc+1))
3956 adv_vert_sg(kc) = ct1_sg(kc)-ct2_sg(kc)-ct3_sg(kc)-ct4_sg(kc)
3957 abs_adv_vert_sg(kc) = abs(adv_vert_sg(kc))
3963 temp_c_help(kc) = 0.5_dp*(
temp_c(kc,j,i)+
temp_c(kc+1,j,i))
3967 ci2(kc) = ai2(kc)/
h_c(j,i)
3971 dtt_dxi = 2.0_dp*dtt_2dxi
3972 dtt_deta = 2.0_dp*dtt_2deta
3979 lgs_a2(kr) = -1.0_dp
3987 lgs_a1(kr) = 1.0_dp + 2.0_dp*ctr1
3989 lgs_b(kr) =
temp_r(kr,j,i)
3998 lgs_a2(kr) = -1.0_dp
3999 lgs_b(kr) = 2.0_dp*clb1
4011 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, krmax)
4030 lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
4032 lgs_a1(kc) = 1.0_dp+ct5(kc)*(ct6(kc)+ct6(kc-1))
4033 lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
4039 = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
4043 +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
4044 -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
4045 +ct5(kc)*(ct6(kc)+ct6(kc-1))
4047 = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) ) &
4052 adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
4055 = -max(adv_vert_help, 0.0_dp) &
4059 +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp) &
4060 +ct5(kc)*(ct6(kc)+ct6(kc-1))
4062 = min(adv_vert_help, 0.0_dp) &
4069 lgs_b(kc) =
temp_c(kc,j,i) + ct7(kc) &
4071 ( (
vx_c(kc,j,i)-abs(
vx_c(kc,j,i))) &
4074 +(
vx_c(kc,j,i-1)+abs(
vx_c(kc,j,i-1))) &
4078 ( (
vy_c(kc,j,i)-abs(
vy_c(kc,j,i))) &
4081 +(
vy_c(kc,j-1,i)+abs(
vy_c(kc,j-1,i))) &
4087 vx_c_help = 0.5_dp*(
vx_c(kc,j,i)+
vx_c(kc,j,i-1))
4088 vy_c_help = 0.5_dp*(
vy_c(kc,j,i)+
vy_c(kc,j-1,i))
4090 lgs_b(kc) =
temp_c(kc,j,i) + ct7(kc) &
4092 ( min(vx_c_help, 0.0_dp) &
4095 +max(vx_c_help, 0.0_dp) &
4099 ( min(vy_c_help, 0.0_dp) &
4102 +max(vy_c_help, 0.0_dp) &
4117 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
4139 lgs_a1(kc) = 1.0_dp - min(adv_vert_sg(kc), 0.0_dp)
4140 lgs_a2(kc) = min(adv_vert_sg(kc), 0.0_dp)
4144 lgs_b(kc) =
age_c(kc,j,i) + dtime_temp &
4146 ( (
vx_c(kc,j,i)-abs(
vx_c(kc,j,i))) &
4149 +(
vx_c(kc,j,i-1)+abs(
vx_c(kc,j,i-1))) &
4153 ( (
vy_c(kc,j,i)-abs(
vy_c(kc,j,i))) &
4156 +(
vy_c(kc,j-1,i)+abs(
vy_c(kc,j-1,i))) &
4162 vx_c_help = 0.5_dp*(
vx_c(kc,j,i)+
vx_c(kc,j,i-1))
4163 vy_c_help = 0.5_dp*(
vy_c(kc,j,i)+
vy_c(kc,j-1,i))
4165 lgs_b(kc) =
age_c(kc,j,i) + dtime_temp &
4167 ( min(vx_c_help, 0.0_dp) &
4170 +max(vx_c_help, 0.0_dp) &
4174 ( min(vy_c_help, 0.0_dp) &
4177 +max(vy_c_help, 0.0_dp) &
4187 lgs_a0(kc) = -0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
4189 lgs_a1(kc) = 1.0_dp+ci1(kc)*(ci2(kc)+ci2(kc-1))
4190 lgs_a2(kc) = 0.5_dp*(ct1(kc)-ct2(kc)-ct3(kc)-ct4(kc)) &
4195 lgs_a0(kc) = -0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1))
4196 lgs_a1(kc) = 1.0_dp &
4197 +0.5_dp*(adv_vert_sg(kc-1)+abs_adv_vert_sg(kc-1)) &
4198 -0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
4199 lgs_a2(kc) = 0.5_dp*(adv_vert_sg(kc) -abs_adv_vert_sg(kc) )
4203 adv_vert_help = 0.5_dp*(adv_vert_sg(kc)+adv_vert_sg(kc-1))
4205 lgs_a0(kc) = -max(adv_vert_help, 0.0_dp)
4206 lgs_a1(kc) = 1.0_dp &
4207 +max(adv_vert_help, 0.0_dp)-min(adv_vert_help, 0.0_dp)
4208 lgs_a2(kc) = min(adv_vert_help, 0.0_dp)
4214 lgs_b(kc) =
age_c(kc,j,i) + dtime_temp &
4216 ( (
vx_c(kc,j,i)-abs(
vx_c(kc,j,i))) &
4219 +(
vx_c(kc,j,i-1)+abs(
vx_c(kc,j,i-1))) &
4223 ( (
vy_c(kc,j,i)-abs(
vy_c(kc,j,i))) &
4226 +(
vy_c(kc,j-1,i)+abs(
vy_c(kc,j-1,i))) &
4232 vx_c_help = 0.5_dp*(
vx_c(kc,j,i)+
vx_c(kc,j,i-1))
4233 vy_c_help = 0.5_dp*(
vy_c(kc,j,i)+
vy_c(kc,j-1,i))
4235 lgs_b(kc) =
age_c(kc,j,i) + dtime_temp &
4237 ( min(vx_c_help, 0.0_dp) &
4240 +max(vx_c_help, 0.0_dp) &
4244 ( min(vy_c_help, 0.0_dp) &
4247 +max(vy_c_help, 0.0_dp) &
4256 if (
as_perp(j,i) >= zero)
then 4261 lgs_a0(kc) = -max(adv_vert_sg(kc-1), 0.0_dp)
4262 lgs_a1(kc) = 1.0_dp + max(adv_vert_sg(kc-1), 0.0_dp)
4267 lgs_b(kc) =
age_c(kc,j,i) + dtime_temp &
4269 ( (
vx_c(kc,j,i)-abs(
vx_c(kc,j,i))) &
4272 +(
vx_c(kc,j,i-1)+abs(
vx_c(kc,j,i-1))) &
4276 ( (
vy_c(kc,j,i)-abs(
vy_c(kc,j,i))) &
4279 +(
vy_c(kc,j-1,i)+abs(
vy_c(kc,j-1,i))) &
4285 vx_c_help = 0.5_dp*(
vx_c(kc,j,i)+
vx_c(kc,j,i-1))
4286 vy_c_help = 0.5_dp*(
vy_c(kc,j,i)+
vy_c(kc,j-1,i))
4288 lgs_b(kc) =
age_c(kc,j,i) + dtime_temp &
4290 ( min(vx_c_help, 0.0_dp) &
4293 +max(vx_c_help, 0.0_dp) &
4297 ( min(vy_c_help, 0.0_dp) &
4300 +max(vy_c_help, 0.0_dp) &
4310 call tri_sle(lgs_a0, lgs_a1, lgs_a2, lgs_x, lgs_b, kcmax)
4319 if (
age_c_neu(kc,j,i) < (age_min*year_sec)) &
4321 if (
age_c_neu(kc,j,i) > (age_max*year_sec)) &
4332 end subroutine calc_temp_ssa
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) vx_c
vx_c(kc,j,i): Velocity in x-direction in the upper (kc) ice domain (at (i+1/2,j,kc)) ...
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) temp_c
temp_c(kc,j,i): Temperature in the upper (kc) ice domain
real(dp) kappa_r
KAPPA_R: Heat conductivity of the lithosphere.
real(dp) delta_tm_sw
DELTA_TM_SW: Melting point depression of sea water due to its average salinity.
integer(i1b), dimension(0:jmax, 0:imax) n_cts
n_cts(j,i): Mask for thermal conditions. -1: cold ice base, 0: temperate ice base with cold ice above...
real(dp), dimension(0:krmax, 0:jmax, 0:imax) temp_r_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) de_c
de_c(kc,j,i): Full effective strain rate in the upper (kc) ice domain
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) age_c
age_c(kc,j,i): Age in the upper (kc) ice domain
real(dp), dimension(0:jmax, 0:imax) insq_g22_sgy
insq_g22_sgy(j,i): Inverse square root of g22, at (i,j+1/2)
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) de_t
de_t(kt,j,i): Full effective strain rate in the lower (kt) ice domain
real(dp), dimension(0:kcmax) eaz_c
eaz_c(kc): Abbreviation for exp(aa*zeta(kc))
real(dp), dimension(0:jmax, 0:imax) q_geo
q_geo(j,i): Geothermal heat flux
real(dp), dimension(0:jmax, 0:imax) vx_b_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) age_t
age_t(kt,j,i): Age in the lower (kt) ice domain
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) enh_c
enh_c(kc,j,i): Flow enhancement factor in the upper (kc) ice domain
real(dp), dimension(0:krmax, 0:jmax, 0:imax) temp_r
temp_r(kr,j,i): Temperature in the bedrock
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) age_t_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp) h_r
H_R: Thickness of the modelled lithosphere layer.
subroutine error(error_message)
Main routine of error_m: Writing of error messages and stopping execution.
Material properties of ice: Rate factor, heat conductivity, specific heat (heat capacity), creep function, viscosity.
subroutine calc_temp_r(atr1, alb1, i, j)
Computation of temperature and age for an ice-free column.
real(dp), parameter eps
eps: Small number
real(dp) omega_max
OMEGA_MAX: Threshold value for the water content of temperate ice.
real(dp), dimension(0:jmax, 0:imax) insq_g11_sgx
insq_g11_sgx(j,i): Inverse square root of g11, at (i+1/2,j)
Several mathematical tools used by SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) insq_g22_g
insq_g22_g(j,i): Inverse square root of g22 on grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) dh_t_dtau
dH_t_dtau(j,i): Derivative of H_t by tau (time)
real(dp), dimension(0:jmax, 0:imax) dh_c_dtau
dH_c_dtau(j,i): Derivative of H_c by tau (time)
real(dp) function, public kappa_val(temp_val)
Heat conductivity of ice: Linear interpolation of tabulated values in KAPPA(.).
real(dp), dimension(0:jmax, 0:imax) dh_c_dxi_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
subroutine, public calc_temp_cold(dxi, deta, dzeta_c, dzeta_t, dzeta_r, dtime_temp)
Computation of temperature and age in cold-ice mode.
subroutine calc_temp2(at1, at2_1, at2_2, at3_1, at3_2, at4_1, at4_2, at5, at6, at7, atr1, alb1, ai1, ai2, dtime_temp, dtt_2dxi, dtt_2deta, i, j)
Computation of temperature and age for an ice column with a temperate base overlain by cold ice...
real(dp), dimension(0:jmax, 0:imax) as_perp
as_perp(j,i): Accumulation-ablation function at the ice surface (SMB)
real(dp), dimension(0:jmax, 0:imax) area
area(j,i): Area of grid cell associated with grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) dzs_deta_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) temp_t_m
temp_t_m(kt,j,i): Melting temperature in the lower (kt) ice domain
real(dp), dimension(0:jmax, 0:imax) de_ssa
de_ssa(j,i): Effective strain rate of the SSA, at (i,j)
real(dp), dimension(0:jmax, 0:imax) h_t_neu
(.)_neu: New value of quantity (.) computed during an integration step
subroutine, public calc_temp_const()
Isothermal mode: Setting of the temperature and age to constant values.
subroutine tri_sle(a0, a1, a2, x, b, nrows)
Solution of a system of linear equations Ax=b with tridiagonal matrix A.
real(dp) function, public creep(sigma_val)
Creep response function for ice.
real(dp), dimension(0:jmax, 0:imax) temp_s
temp_s(j,i): Ice surface temperature
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) vz_t
vz_t(kt,j,i): Velocity in z-direction in the lower (kt) ice domain (at (i,j,kt+1/2)) ...
real(dp) g
G: Acceleration due to gravity.
real(dp), dimension(0:jmax, 0:imax) q_tld
Q_tld(j,i): Water drainage rate from the temperate layer.
Declarations of global variables for SICOPOLIS (for the ANT domain).
subroutine calc_temp3(at1, at2_1, at2_2, at3_1, at3_2, at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, ai1, ai2, ai3, dzeta_t, dtime_temp, dtt_2dxi, dtt_2deta, i, j)
Computation of temperature, water content and age for an ice column with a temperate base overlain by...
real(dp) function, public ratefac_c(temp_val, temp_m_val)
Rate factor for cold ice: Linear interpolation of tabulated values in RF(.).
real(dp) ea
ea: Abbreviation for exp(aa)
logical flag_aa_nonzero
flag_aa_nonzero: Flag for the exponential stretch parameter aa. .true.: aa greater than zero (non-equ...
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) age_c_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp) aa
aa: Exponential stretch parameter of the non-equidistant vertical grid in the upper (kc) ice domain ...
real(dp) rho_c_r
RHO_C_R: Density times specific heat of the lithosphere.
Declarations of kind types for SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) zs
zs(j,i): Coordinate z of the surface topography
real(dp), dimension(0:jmax, 0:imax) dh_c_deta_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) enh_t
enh_t(kt,j,i): Flow enhancement factor in the lower (kt) ice domain
real(dp) function, public ratefac_t(omega_val)
Rate factor for temperate ice.
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) temp_c_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp), dimension(0:jmax, 0:imax) h_c
H_c(j,i): Thickness of ice in the upper (kc) domain (thickness of the cold-ice layer for POLY...
Computation of temperature, water content and age.
real(dp), dimension(0:jmax, 0:imax) vy_b_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
integer(i1b), dimension(0:jmax, 0:imax) n_cts_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) omega_t_neu
(.)_neu: New value of quantity (.) computed during an integration step
Writing of error messages and stopping execution.
real(dp), dimension(0:jmax, 0:imax) dzm_dtau
dzm_dtau(j,i): Derivative of zm by tau (time)
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) temp_c_m
temp_c_m(kc,j,i): Melting temperature in the upper (kc) ice domain
real(dp), dimension(0:jmax, 0:imax) dzb_deta_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp) nue
NUE: Water diffusivity in ice.
real(dp) function, public c_val(temp_val)
Specific heat of ice: Linear interpolation of tabulated values in C(.).
real(dp), dimension(0:jmax, 0:imax) p_weert_inv
p_weert_inv(j,i): Inverse of p_weert
logical, dimension(0:jmax, 0:imax) flag_shelfy_stream
flag_shelfy_stream(j,i): Shelfy stream flag on the main grid. .true.: grounded ice, and at least one neighbour on the staggered grid is a shelfy stream point .false.: otherwise
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) vz_c
vz_c(kc,j,i): Velocity in z-direction in the upper (kc) ice domain (at (i,j,kc+1/2)) ...
integer(i4b), dimension(0:jmax, 0:imax) kc_cts_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp), dimension(0:jmax, 0:imax) h_c_neu
(.)_neu: New value of quantity (.) computed during an integration step
real(dp), dimension(0:jmax, 0:imax) zm
zm(j,i): Coordinate z of the bottom of the upper (kc) ice domain = top of the lower (kt) ice domain (...
integer(i1b), dimension(0:jmax, 0:imax) maske
maske(j,i): Ice-land-ocean mask. 0: grounded ice, 1: ice-free land, 2: ocean, 3: floating ice ...
real(dp), dimension(0:kcmax) eaz_c_quotient
eaz_c_quotient(kc): Abbreviation for (eaz_c(kc)-1.0)/(ea-1.0)
real(dp), dimension(0:jmax, 0:imax) dh_t_deta_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) dh_t_dxi_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
subroutine shift_cts_upward(at1, at2_1, at2_2, at3_1, at3_2, at4_1, at4_2, at5, at6, at7, atr1, am1, am2, alb1, aw1, aw2, aw3, aw4, aw5, aw7, aw8, aw9, aqtld, ai1, ai2, ai3, dzeta_t, dtime_temp, dtt_2dxi, dtt_2deta, dtime_temp_inv, i, j)
Upward shifting of the CTS.
real(dp), dimension(0:jmax, 0:imax) dzm_dxi_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) dzs_dxi_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) vy_t
vy_t(kt,j,i): Velocity in y-direction in the lower (kt) ice domain (at (i,j+1/2,kt)) ...
real(dp), dimension(0:jmax, 0:imax) am_perp
am_perp(j,i): Ice volume flux across the z=zm interface
real(dp) beta
BETA: Clausius-Clapeyron gradient of ice.
real(dp), dimension(0:jmax, 0:imax) dzm_deta_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) insq_g11_g
insq_g11_g(j,i): Inverse square root of g11 on grid point (i,j)
real(dp), dimension(0:jmax, 0:imax) c_drag
c_drag(j,i): Auxiliary quantity for the computation of the basal drag
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) vy_c
vy_c(kc,j,i): Velocity in y-direction in the upper (kc) ice domain (at (i,j+1/2,kc)) ...
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) omega_t
omega_t(kt,j,i): Water content in the lower (kt) ice domain
real(dp), dimension(0:ktmax) zeta_t
zeta_t(kt): Sigma coordinate zeta_t of grid point kt
real(dp) function, public viscosity(de_val, temp_val, temp_m_val, omega_val, enh_val, i_flag_cold_temp)
Ice viscosity as a function of the effective strain rate and the temperature (in cold ice) or the wat...
real(dp) l
L: Latent heat of ice.
real(dp), dimension(0:jmax, 0:imax) zm_neu
(.)_neu: New value of quantity (.) computed during an integration step
character(len=256) errormsg
errormsg: Error-message string
real(dp) rho
RHO: Density of ice.
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) sigma_c
sigma_c(kc,j,i): Effective stress in the upper (kc) ice domain
real(dp), dimension(0:jmax, 0:imax) zb
zb(j,i): Coordinate z of the ice base
real(dp), dimension(0:kcmax, 0:jmax, 0:imax) omega_c_neu
(.)_neu: New value of quantity (.) computed during an integration step
subroutine, public calc_temp_poly(dxi, deta, dzeta_c, dzeta_t, dzeta_r, dtime_temp)
Computation of temperature, water content and age in polythermal mode.
real(dp), dimension(0:jmax, 0:imax) am_perp_st
am_perp_st(j,i): Steady-state part of am_perp (without contribution of dzm_dtau)
real(dp), dimension(0:jmax, 0:imax) dzb_dtau
dzb_dtau(j,i): Derivative of zb by tau (time)
real(dp), dimension(0:kcmax) zeta_c
zeta_c(kc): Sigma coordinate zeta_c of grid point kc
real(dp), dimension(0:jmax, 0:imax) h_t
H_t(j,i): Thickness of ice in the lower (kt) domain (thickness of the temperate layer for POLY...
Declarations of global variables for SICOPOLIS.
real(dp), dimension(0:jmax, 0:imax) dzs_dtau
dzs_dtau(j,i): Derivative of zs by tau (time)
real(dp), dimension(0:ktmax, 0:jmax, 0:imax) vx_t
vx_t(kt,j,i): Velocity in x-direction in the lower (kt) ice domain (at (i+1/2,j,kt)) ...
real(dp), dimension(0:jmax, 0:imax) dzb_dxi_g
(.)_g(j,i): Staggered-grid quantity (.) interpolated to grid point (i,j)