1 | MODULE soil_evolution_mod |
---|
2 | |
---|
3 | IMPLICIT NONE |
---|
4 | |
---|
5 | LOGICAL soil_pem ! True by default, to run with the subsurface physic. Read in pem.def |
---|
6 | |
---|
7 | CONTAINS |
---|
8 | |
---|
9 | subroutine soil_pem_CN(ngrid,nsoil, & |
---|
10 | therm_i, timestep,tsurf,tsurf_prev,tsoil) |
---|
11 | |
---|
12 | |
---|
13 | use comsoil_h_PEM, only: & fluxgeo,layer_PEM |
---|
14 | use comsoil_h,only: volcapa |
---|
15 | implicit none |
---|
16 | |
---|
17 | !----------------------------------------------------------------------- |
---|
18 | ! Author: LL |
---|
19 | ! Purpose: Compute soil temperature using an implict 1st order scheme |
---|
20 | ! |
---|
21 | ! Note: depths of layers and mid-layers, soil thermal inertia and |
---|
22 | ! heat capacity are commons in comsoil_PEM.h |
---|
23 | ! A convergence loop is added until equilibrium |
---|
24 | !----------------------------------------------------------------------- |
---|
25 | |
---|
26 | #include "dimensions.h" |
---|
27 | |
---|
28 | !----------------------------------------------------------------------- |
---|
29 | ! arguments |
---|
30 | ! --------- |
---|
31 | ! inputs: |
---|
32 | integer,intent(in) :: ngrid ! number of (horizontal) grid-points |
---|
33 | integer,intent(in) :: nsoil ! number of soil layers |
---|
34 | real,intent(in) :: therm_i(ngrid,nsoil) ! thermal inertia |
---|
35 | real,intent(in) :: timestep ! time step |
---|
36 | real,intent(in) :: tsurf(ngrid) ! surface temperature |
---|
37 | real,intent(in) :: tsurf_prev(ngrid) ! surface temperature at previous timestep |
---|
38 | |
---|
39 | ! outputs: |
---|
40 | real,intent(inout) :: tsoil(ngrid,nsoil) ! soil (mid-layer) temperature |
---|
41 | |
---|
42 | ! local variables: |
---|
43 | integer ig,ik,isoil |
---|
44 | REAL :: ice_inertia = 2120. |
---|
45 | real :: alph_PEM(nsoil) |
---|
46 | real :: beta_PEM(nsoil) |
---|
47 | real :: rhoc(nsoil) |
---|
48 | real :: volcapa_ice = 1.43e7 |
---|
49 | |
---|
50 | real :: k_soil(nsoil) |
---|
51 | real :: d1(nsoil),d2(nsoil),d3(nsoil),re(nsoil) |
---|
52 | real :: Tcol(nsoil) |
---|
53 | |
---|
54 | do ig = 1,ngrid |
---|
55 | |
---|
56 | Tcol(:) = tsoil(ig,:) |
---|
57 | !0. Build thermal properties of the ground |
---|
58 | |
---|
59 | do isoil = 1,nsoil |
---|
60 | if(abs(ice_inertia-therm_i(ig,isoil)).lt.1e-6) then |
---|
61 | rhoc(isoil) = volcapa_ice |
---|
62 | else |
---|
63 | rhoc(isoil) = volcapa |
---|
64 | endif |
---|
65 | k_soil(isoil) = therm_i(ig,isoil)*therm_i(ig,isoil)/rhoc(isoil) |
---|
66 | enddo |
---|
67 | |
---|
68 | !1. Build Crank-Nicholson scheme for heat equation, following Schorgofer derivations |
---|
69 | !1.a all the points which are not impacted by boundary condition |
---|
70 | do isoil = 2,nsoil-1 |
---|
71 | |
---|
72 | alph_pem(isoil) = timestep/((rhoc(isoil) + rhoc(isoil+1))/2.)*k_soil(isoil+1)/((layer_PEM(isoil+1) - layer_PEM(isoil))*(layer_PEM(isoil+1) - layer_PEM(isoil-1))) |
---|
73 | beta_pem(isoil) = timestep/((rhoc(isoil) + rhoc(isoil+1))/2.)*k_soil(isoil)/((layer_PEM(isoil) - layer_PEM(isoil-1))*(layer_PEM(isoil+1) - layer_PEM(isoil-1))) |
---|
74 | enddo |
---|
75 | |
---|
76 | !1.b At the surface: T(z=0) = Tsurf. |
---|
77 | |
---|
78 | alph_pem(1) = k_soil(2)*timestep/(rhoc(1)*layer_PEM(2)*(layer_PEM(2)-layer_PEM(1))) |
---|
79 | beta_pem(1) = k_soil(1)*timestep/(rhoc(1)*layer_PEM(2)*layer_PEM(1)) |
---|
80 | |
---|
81 | !1.c At the bottom: dT/dz = -Fgeo/k. |
---|
82 | beta_pem(nsoil) = timestep*k_soil(nsoil)/(2.*rhoc(nsoil)*(layer_PEM(isoil) - layer_PEM(isoil-1))**2) |
---|
83 | |
---|
84 | !1 .d set the tridiagonal with diagonals d1,d2,d3, and right element re: |
---|
85 | |
---|
86 | |
---|
87 | d1(:) = -beta_pem(:) |
---|
88 | d2(:) = 1. + alph_pem(:) + beta_pem(:) |
---|
89 | d3(:) = -alph_pem(:) |
---|
90 | d2(nsoil) = 1. + beta_pem(nsoil) |
---|
91 | |
---|
92 | re(1)= alph_pem(1)*Tcol(2) + (1.-alph_pem(1)-beta_pem(1))*Tcol(1) + beta_pem(1)*(tsurf(ig)+tsurf_prev(ig)) |
---|
93 | do isoil=2,nsoil-1 |
---|
94 | re(isoil) = beta_pem(isoil)*Tcol(isoil-1) + (1.-alph_pem(isoil)-beta_pem(isoil))*Tcol(isoil) + alph_pem(isoil)*Tcol(isoil+1) |
---|
95 | enddo |
---|
96 | re(nsoil) = beta_pem(nsoil)*Tcol(nsoil-1) + (1.-beta_pem(nsoil))*Tcol(nsoil) + & |
---|
97 | timestep/rhoc(nsoil)*fluxgeo/(layer_PEM(nsoil)-layer_PEM(nsoil-1)) |
---|
98 | |
---|
99 | !2. Update by tridiagonal inversion |
---|
100 | call tridag(d1,d2,d3,re,Tcol,nsoil) |
---|
101 | tsoil(ig,:) = Tcol(:) |
---|
102 | enddo |
---|
103 | |
---|
104 | |
---|
105 | END SUBROUTINE soil_pem_CN |
---|
106 | |
---|
107 | !==================================== |
---|
108 | |
---|
109 | SUBROUTINE TRIDAG(A,B,C,R,U,N) |
---|
110 | ! From Numerical Recipes in Fortran 77 |
---|
111 | INTEGER, INTENT(IN):: N |
---|
112 | INTEGER ,PARAMETER:: NMAX = 1000 |
---|
113 | |
---|
114 | REAL, INTENT(IN) :: A(N),B(N),C(N),R(N) |
---|
115 | REAL,INTENT(INOUT) :: U(N) |
---|
116 | |
---|
117 | INTEGER j |
---|
118 | REAL bet,gam(NMAX) |
---|
119 | |
---|
120 | IF(B(1).EQ.0.)PAUSE |
---|
121 | BET=B(1) |
---|
122 | U(1)=R(1)/BET |
---|
123 | DO 11 J=2,N |
---|
124 | GAM(J)=C(J-1)/BET |
---|
125 | BET=B(J)-A(J)*GAM(J) |
---|
126 | IF(BET.EQ.0.)PAUSE |
---|
127 | U(J)=(R(J)-A(J)*U(J-1))/BET |
---|
128 | 11 CONTINUE |
---|
129 | DO 12 J=N-1,1,-1 |
---|
130 | U(J)=U(J)-GAM(J+1)*U(J+1) |
---|
131 | 12 CONTINUE |
---|
132 | RETURN |
---|
133 | END SUBROUTINE TRIDAG |
---|
134 | |
---|
135 | |
---|
136 | END MODULE soil_evolution_mod |
---|
137 | |
---|
138 | |
---|
139 | |
---|