/[lmdze]/trunk/Sources/phylmd/Interface_surf/soil.f
ViewVC logotype

Contents of /trunk/Sources/phylmd/Interface_surf/soil.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 157 - (show annotations)
Mon Jul 20 16:01:49 2015 UTC (8 years, 10 months ago) by guez
File size: 7848 byte(s)
Just encapsulated SUBROUTINE vlsplt in a module and cleaned it.

In procedure vlx, local variables dxqu and adxqu only need indices
iip2:ip1jm. Otherwise, just cleaned vlx.

Procedures dynredem0 and dynredem1 no longer have argument fichnom,
they just operate on a file named "restart.nc". The programming
guideline here is that gcm should not be more complex than it needs by
itself, other programs (ce0l etc.) just have to adapt to gcm. So ce0l
now creates files "restart.nc" and "restartphy.nc".

In order to facilitate decentralizing the writing of "restartphy.nc",
created a procedure phyredem0 out of phyredem. phyredem0 creates the
NetCDF header of "restartphy.nc" while phyredem writes the NetCDF
variables. As the global attribute itau_phy needs to be filled in
phyredem0, at the beginnig of the run, we must compute its value
instead of just using itap. So we have a dummy argument lmt_pas of
phyredem0. Also, the ncid of "startphy.nc" is upgraded from local
variable of phyetat0 to dummy argument. phyetat0 no longer closes
"startphy.nc".

Following the same decentralizing objective, the ncid of "restart.nc"
is upgraded from local variable of dynredem0 to module variable of
dynredem0_m. "restart.nc" is not closed at the end of dynredem0 nor
opened at the beginning of dynredem1.

In procedure etat0, instead of creating many vectors of size klon
which will be filled with zeroes, just create one array null_array.

In procedure phytrac, instead of writing trs(: 1) to a text file,
write it to "restartphy.nc" (following LMDZ). This is better because
now trs(: 1) is next to its coordinates. We can write to
"restartphy.nc" from phytrac directly, and not add trs(: 1) to the
long list of variables in physiq, thanks to the decentralizing of
"restartphy.nc".

In procedure phyetat0, we no longer write to standard output the
minimum and maximum values of read arrays. It is ok to check input and
abort on invalid values but just printing statistics on input seems too
much useless computation and out of place clutter.

1 module soil_m
2
3 IMPLICIT NONE
4
5 contains
6
7 SUBROUTINE soil(ptimestep, indice, knon, snow, ptsrf, ptsoil, pcapcal, &
8 pfluxgrd)
9
10 ! From LMDZ4/libf/phylmd/soil.F, version 1.1.1.1 2004/05/19
11
12 USE dimens_m
13 USE indicesol
14 USE dimphy
15 USE dimsoil
16 USE suphec_m
17
18 ! =======================================================================
19
20 ! Auteur: Frederic Hourdin 30/01/92
21 ! -------
22
23 ! objet: computation of : the soil temperature evolution
24 ! ------ the surfacic heat capacity "Capcal"
25 ! the surface conduction flux pcapcal
26
27
28 ! Method: implicit time integration
29 ! -------
30 ! Consecutive ground temperatures are related by:
31 ! T(k+1) = C(k) + D(k)*T(k) (1)
32 ! the coefficients C and D are computed at the t-dt time-step.
33 ! Routine structure:
34 ! 1)new temperatures are computed using (1)
35 ! 2)C and D coefficients are computed from the new temperature
36 ! profile for the t+dt time-step
37 ! 3)the coefficients A and B are computed where the diffusive
38 ! fluxes at the t+dt time-step is given by
39 ! Fdiff = A + B Ts(t+dt)
40 ! or Fdiff = F0 + Capcal (Ts(t+dt)-Ts(t))/dt
41 ! with F0 = A + B (Ts(t))
42 ! Capcal = B*dt
43
44 ! Interface:
45 ! ----------
46
47 ! Arguments:
48 ! ----------
49 ! ptimestep physical timestep (s)
50 ! indice sub-surface index
51 ! snow(klon,nbsrf) snow
52 ! ptsrf(knon) surface temperature at time-step t (K)
53 ! ptsoil(klon,nsoilmx) temperature inside the ground (K)
54 ! pcapcal(klon) surfacic specific heat (W*m-2*s*K-1)
55 ! pfluxgrd(klon) surface diffusive flux from ground (Wm-2)
56
57 ! =======================================================================
58 ! declarations:
59 ! -------------
60
61
62 ! -----------------------------------------------------------------------
63 ! arguments
64 ! ---------
65
66 REAL ptimestep
67 INTEGER indice, knon
68 REAL ptsrf(knon), ptsoil(klon, nsoilmx), snow(klon)
69 REAL pcapcal(klon), pfluxgrd(klon)
70
71 ! -----------------------------------------------------------------------
72 ! local arrays
73 ! ------------
74
75 INTEGER ig, jk
76 ! $$$ REAL zdz2(nsoilmx),z1(klon)
77 REAL zdz2(nsoilmx), z1(klon, nbsrf)
78 REAL min_period, dalph_soil
79 REAL ztherm_i(klon)
80
81 ! local saved variables:
82 ! ----------------------
83 REAL dz1(nsoilmx), dz2(nsoilmx)
84 ! $$$ REAL zc(klon,nsoilmx),zd(klon,nsoilmx)
85 REAL zc(klon, nsoilmx, nbsrf), zd(klon, nsoilmx, nbsrf)
86 REAL lambda
87 SAVE dz1, dz2, zc, zd, lambda
88 LOGICAL firstsurf(nbsrf)
89 SAVE firstsurf
90 REAL isol, isno, iice
91 SAVE isol, isno, iice
92
93 DATA firstsurf/.TRUE., .TRUE., .TRUE., .TRUE./
94
95 DATA isol, isno, iice/2000., 2000., 2000./
96
97 ! -----------------------------------------------------------------------
98 ! Depthts:
99 ! --------
100
101 REAL rk, fz1, rk1, rk2
102
103 pfluxgrd(:) = 0.
104 ! calcul de l'inertie thermique a partir de la variable rnat.
105 ! on initialise a iice meme au-dessus d'un point de mer au cas
106 ! ou le point de mer devienne point de glace au pas suivant
107 ! on corrige si on a un point de terre avec ou sans glace
108
109 IF (indice==is_sic) THEN
110 DO ig = 1, knon
111 ztherm_i(ig) = iice
112 IF (snow(ig)>0.0) ztherm_i(ig) = isno
113 END DO
114 ELSE IF (indice==is_lic) THEN
115 DO ig = 1, knon
116 ztherm_i(ig) = iice
117 IF (snow(ig)>0.0) ztherm_i(ig) = isno
118 END DO
119 ELSE IF (indice==is_ter) THEN
120 DO ig = 1, knon
121 ztherm_i(ig) = isol
122 IF (snow(ig)>0.0) ztherm_i(ig) = isno
123 END DO
124 ELSE IF (indice==is_oce) THEN
125 DO ig = 1, knon
126 ztherm_i(ig) = iice
127 END DO
128 ELSE
129 PRINT *, 'valeur d indice non prevue', indice
130 STOP 1
131 END IF
132
133 IF (firstsurf(indice)) THEN
134
135 ! -----------------------------------------------------------------------
136 ! ground levels
137 ! grnd=z/l where l is the skin depth of the diurnal cycle:
138 ! --------------------------------------------------------
139
140 min_period = 1800. ! en secondes
141 dalph_soil = 2. ! rapport entre les epaisseurs de 2 couches succ.
142
143 OPEN (99, FILE='soil.def', STATUS='old', FORM='formatted', ERR=9999)
144 READ (99, *) min_period
145 READ (99, *) dalph_soil
146 PRINT *, 'Discretization for the soil model'
147 PRINT *, 'First level e-folding depth', min_period, ' dalph', &
148 dalph_soil
149 CLOSE (99)
150 9999 CONTINUE
151
152 ! la premiere couche represente un dixieme de cycle diurne
153 fz1 = sqrt(min_period/3.14)
154
155 DO jk = 1, nsoilmx
156 rk1 = jk
157 rk2 = jk - 1
158 dz2(jk) = fz(rk1) - fz(rk2)
159 END DO
160 DO jk = 1, nsoilmx - 1
161 rk1 = jk + .5
162 rk2 = jk - .5
163 dz1(jk) = 1./(fz(rk1)-fz(rk2))
164 END DO
165 lambda = fz(.5)*dz1(1)
166 PRINT *, 'full layers, intermediate layers (seconds)'
167 DO jk = 1, nsoilmx
168 rk = jk
169 rk1 = jk + .5
170 rk2 = jk - .5
171 PRINT *, 'fz=', fz(rk1)*fz(rk2)*3.14, fz(rk)*fz(rk)*3.14
172 END DO
173 ! PB
174 firstsurf(indice) = .FALSE.
175
176 ! Initialisations:
177 ! ----------------
178
179 ELSE
180 ! -----------------------------------------------------------------------
181 ! Computation of the soil temperatures using the Cgrd and Dgrd
182 ! coefficient computed at the previous time-step:
183 ! -----------------------------------------------
184
185 ! surface temperature
186 DO ig = 1, knon
187 ptsoil(ig, 1) = (lambda*zc(ig,1,indice)+ptsrf(ig))/(lambda*(1.-zd(ig,1, &
188 indice))+1.)
189 END DO
190
191 ! other temperatures
192 DO jk = 1, nsoilmx - 1
193 DO ig = 1, knon
194 ptsoil(ig, jk+1) = zc(ig, jk, indice) + zd(ig, jk, indice)*ptsoil(ig, &
195 jk)
196 END DO
197 END DO
198
199 END IF
200 ! -----------------------------------------------------------------------
201 ! Computation of the Cgrd and Dgrd coefficient for the next step:
202 ! ---------------------------------------------------------------
203
204 ! $$$ PB ajout pour cas glace de mer
205 IF (indice==is_sic) THEN
206 DO ig = 1, knon
207 ptsoil(ig, nsoilmx) = rtt - 1.8
208 END DO
209 END IF
210
211 DO jk = 1, nsoilmx
212 zdz2(jk) = dz2(jk)/ptimestep
213 END DO
214
215 DO ig = 1, knon
216 z1(ig, indice) = zdz2(nsoilmx) + dz1(nsoilmx-1)
217 zc(ig, nsoilmx-1, indice) = zdz2(nsoilmx)*ptsoil(ig, nsoilmx)/ &
218 z1(ig, indice)
219 zd(ig, nsoilmx-1, indice) = dz1(nsoilmx-1)/z1(ig, indice)
220 END DO
221
222 DO jk = nsoilmx - 1, 2, -1
223 DO ig = 1, knon
224 z1(ig, indice) = 1./(zdz2(jk)+dz1(jk-1)+dz1(jk)*(1.-zd(ig,jk,indice)))
225 zc(ig, jk-1, indice) = (ptsoil(ig,jk)*zdz2(jk)+dz1(jk)*zc(ig,jk,indice) &
226 )*z1(ig, indice)
227 zd(ig, jk-1, indice) = dz1(jk-1)*z1(ig, indice)
228 END DO
229 END DO
230
231 ! -----------------------------------------------------------------------
232 ! computation of the surface diffusive flux from ground and
233 ! calorific capacity of the ground:
234 ! ---------------------------------
235
236 DO ig = 1, knon
237 pfluxgrd(ig) = ztherm_i(ig)*dz1(1)*(zc(ig,1,indice)+(zd(ig,1, &
238 indice)-1.)*ptsoil(ig,1))
239 pcapcal(ig) = ztherm_i(ig)*(dz2(1)+ptimestep*(1.-zd(ig,1,indice))*dz1(1))
240 z1(ig, indice) = lambda*(1.-zd(ig,1,indice)) + 1.
241 pcapcal(ig) = pcapcal(ig)/z1(ig, indice)
242 pfluxgrd(ig) = pfluxgrd(ig) + pcapcal(ig)*(ptsoil(ig,1)*z1(ig,indice)- &
243 lambda*zc(ig,1,indice)-ptsrf(ig))/ptimestep
244 END DO
245
246 contains
247
248 real function fz(rk)
249 real rk
250 fz = fz1*(dalph_soil**rk-1.)/(dalph_soil-1.)
251 end function fz
252
253 END SUBROUTINE soil
254
255 end module soil_m

  ViewVC Help
Powered by ViewVC 1.1.21