/[lmdze]/trunk/Sources/dyn3d/etat0.f
ViewVC logotype

Contents of /trunk/Sources/dyn3d/etat0.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: 11798 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 etat0_mod
2
3 use indicesol, only: nbsrf
4 use dimphy, only: klon
5
6 IMPLICIT NONE
7
8 REAL pctsrf(klon, nbsrf)
9 ! ("pctsrf(i, :)" is the composition of the surface at horizontal
10 ! position "i")
11
12 private nbsrf, klon
13
14 contains
15
16 SUBROUTINE etat0(phis)
17
18 ! From "etat0_netcdf.F", version 1.3, 2005/05/25 13:10:09
19
20 use caldyn0_m, only: caldyn0
21 use comconst, only: cpp, kappa, iniconst
22 use comgeom, only: aire_2d, apoln, apols, cu_2d, cv_2d, inigeom
23 use conf_gcm_m, only: nday, day_step, iphysiq
24 use dimens_m, only: iim, jjm, llm, nqmx
25 use dimphy, only: zmasq
26 use dimsoil, only: nsoilmx
27 use disvert_m, only: ap, bp, preff, pa, disvert
28 use dynetat0_m, only: day_ref, annee_ref, xprimp025, xprimm025, rlatu1, &
29 rlatu2, rlatu, rlatv, yprimu1, yprimu2, rlonu, rlonv, xprimu, xprimv
30 use dynredem0_m, only: dynredem0
31 use dynredem1_m, only: dynredem1
32 use exner_hyb_m, only: exner_hyb
33 use fxhyp_m, only: fxhyp
34 use fyhyp_m, only: fyhyp
35 use geopot_m, only: geopot
36 use grid_atob, only: grille_m
37 use grid_change, only: init_dyn_phy, dyn_phy
38 use indicesol, only: is_oce, is_sic, is_ter, is_lic, epsfra
39 use iniadvtrac_m, only: iniadvtrac
40 use inifilr_m, only: inifilr
41 use massdair_m, only: massdair
42 use netcdf, only: nf90_nowrite
43 use netcdf95, only: nf95_close, nf95_get_var, nf95_gw_var, nf95_put_var, &
44 nf95_inq_varid, nf95_open
45 use nr_util, only: pi, assert
46 use paramet_m, only: ip1jm, ip1jmp1
47 use phyetat0_m, only: rlat, rlon
48 use phyredem0_m, only: phyredem0, ncid_restartphy
49 use phyredem_m, only: phyredem
50 use q_sat_m, only: q_sat
51 use regr_lat_time_coefoz_m, only: regr_lat_time_coefoz
52 use regr_pr_o3_m, only: regr_pr_o3
53 use startdyn, only: start_init_dyn
54 USE start_init_orog_m, only: start_init_orog, mask
55 use start_init_phys_m, only: start_init_phys
56 use start_inter_3d_m, only: start_inter_3d
57 use temps, only: itau_phy
58 use test_disvert_m, only: test_disvert
59 use unit_nml_m, only: unit_nml
60
61 REAL, intent(out):: phis(:, :) ! (iim + 1, jjm + 1)
62 ! surface geopotential, in m2 s-2
63
64 ! Local:
65
66 REAL, dimension(iim + 1, jjm + 1, llm):: ucov, t3d, teta
67 REAL vcov(iim + 1, jjm, llm)
68
69 REAL q(iim + 1, jjm + 1, llm, nqmx)
70 ! (mass fractions of trace species
71 ! "q(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)"
72 ! and pressure level "pls(i, j, l)".)
73
74 real qsat(iim + 1, jjm + 1, llm) ! mass fraction of saturating water vapor
75 REAL qsolsrf(klon, nbsrf), snsrf(klon, nbsrf)
76 REAL albe(klon, nbsrf), evap(klon, nbsrf)
77 REAL tsoil(klon, nsoilmx, nbsrf)
78 REAL null_array(klon)
79 REAL solsw(klon), sollw(klon)
80 !IM "slab" ocean
81 REAL frugs(klon, nbsrf), agesno(klon, nbsrf)
82 REAL rugmer(klon)
83 real, dimension(iim + 1, jjm + 1):: zmea_2d, zstd_2d, zsig_2d, zgam_2d
84 real, dimension(iim + 1, jjm + 1):: zthe_2d, zpic_2d, zval_2d
85 real, dimension(iim + 1, jjm + 1):: tsol_2d, qsol_2d, ps
86 REAL zmea(klon), zstd(klon)
87 REAL zsig(klon), zgam(klon)
88 REAL zthe(klon)
89 REAL zpic(klon), zval(klon)
90 REAL t_ancien(klon, llm), q_ancien(klon, llm)
91 real clwcon(klon, llm), rnebcon(klon, llm), ratqs(klon, llm)
92
93 ! D\'eclarations pour lecture glace de mer :
94 INTEGER iml_lic, jml_lic
95 INTEGER ncid, varid
96 REAL, pointer:: dlon_lic(:), dlat_lic(:)
97 REAL, ALLOCATABLE:: fraclic(:, :) ! fraction land ice
98 REAL flic_tmp(iim + 1, jjm + 1) ! fraction land ice temporary
99
100 INTEGER l, ji
101
102 REAL pk(iim + 1, jjm + 1, llm) ! fonction d'Exner aux milieux des couches
103 real pks(iim + 1, jjm + 1)
104
105 REAL masse(iim + 1, jjm + 1, llm)
106 REAL phi(iim + 1, jjm + 1, llm)
107 REAL pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
108 REAL w(iim + 1, jjm + 1, llm)
109
110 real sig1(klon, llm) ! section adiabatic updraft
111 real w01(klon, llm) ! vertical velocity within adiabatic updraft
112
113 real pls(iim + 1, jjm + 1, llm)
114 ! (pressure at mid-layer of LMDZ grid, in Pa)
115 ! "pls(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)",
116 ! for layer "l")
117
118 REAL p3d(iim + 1, jjm + 1, llm+1) ! pressure at layer interfaces, in Pa
119 ! ("p3d(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)",
120 ! for interface "l")
121
122 namelist /etat0_nml/ day_ref, annee_ref
123
124 !---------------------------------
125
126 print *, "Call sequence information: etat0"
127
128 print *, "Enter namelist 'etat0_nml'."
129 read(unit=*, nml=etat0_nml)
130 write(unit_nml, nml=etat0_nml)
131
132 CALL iniconst
133
134 ! Construct a grid:
135
136 pa = 5e4
137 CALL disvert
138 call test_disvert
139
140 CALL fyhyp(rlatu, rlatv, rlatu2, yprimu2, rlatu1, yprimu1)
141 CALL fxhyp(xprimm025, rlonv, xprimv, rlonu, xprimu, xprimp025)
142
143 rlatu(1) = pi / 2.
144 rlatu(jjm + 1) = -rlatu(1)
145
146 CALL inigeom
147 CALL inifilr
148
149 rlat(1) = 90.
150 rlat(2:klon-1) = pack(spread(rlatu(2:jjm), 1, iim), .true.) * 180. / pi
151 ! (with conversion to degrees)
152 rlat(klon) = - 90.
153
154 rlon(1) = 0.
155 rlon(2:klon-1) = pack(spread(rlonv(:iim), 2, jjm - 1), .true.) * 180. / pi
156 ! (with conversion to degrees)
157 rlon(klon) = 0.
158
159 call start_init_orog(phis, zmea_2d, zstd_2d, zsig_2d, zgam_2d, zthe_2d, &
160 zpic_2d, zval_2d) ! also compute "mask"
161 call init_dyn_phy ! define the mask "dyn_phy" for distinct grid points
162 zmasq = pack(mask, dyn_phy)
163 PRINT *, 'Masque construit'
164
165 call start_init_phys(tsol_2d, qsol_2d)
166 CALL start_init_dyn(tsol_2d, phis, ps)
167
168 ! Compute pressure on intermediate levels:
169 forall(l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * ps
170 CALL exner_hyb(ps, p3d, pks, pk)
171 call assert(MINVAL(pk) /= MAXVAL(pk), '"pk" should not be constant')
172
173 pls = preff * (pk / cpp)**(1. / kappa)
174 PRINT *, "minval(pls) = ", minval(pls)
175 print *, "maxval(pls) = ", maxval(pls)
176
177 call start_inter_3d('U', rlonv, rlatv, pls, ucov)
178 forall (l = 1: llm) ucov(:iim, :, l) = ucov(:iim, :, l) * cu_2d(:iim, :)
179 ucov(iim+1, :, :) = ucov(1, :, :)
180
181 call start_inter_3d('V', rlonu, rlatu(:jjm), pls(:, :jjm, :), vcov)
182 forall (l = 1: llm) vcov(:iim, :, l) = vcov(:iim, :, l) * cv_2d(:iim, :)
183 vcov(iim + 1, :, :) = vcov(1, :, :)
184
185 call start_inter_3d('TEMP', rlonu, rlatv, pls, t3d)
186 PRINT *, 'minval(t3d) = ', minval(t3d)
187 print *, "maxval(t3d) = ", maxval(t3d)
188
189 teta(:iim, :, :) = t3d(:iim, :, :) * cpp / pk(:iim, :, :)
190 teta(iim + 1, :, :) = teta(1, :, :)
191 DO l = 1, llm
192 teta(:, 1, l) = SUM(aire_2d(:, 1) * teta(:, 1, l)) / apoln
193 teta(:, jjm + 1, l) = SUM(aire_2d(:, jjm + 1) * teta(:, jjm + 1, l)) &
194 / apols
195 ENDDO
196
197 ! Calcul de l'humidit\'e \`a saturation :
198 qsat = q_sat(t3d, pls)
199 PRINT *, "minval(qsat) = ", minval(qsat)
200 print *, "maxval(qsat) = ", maxval(qsat)
201 IF (MINVAL(qsat) == MAXVAL(qsat)) stop '"qsat" should not be constant'
202
203 ! Water vapor:
204 call start_inter_3d('R', rlonu, rlatv, pls, q(:, :, :, 1))
205 q(:, :, :, 1) = 0.01 * q(:, :, :, 1) * qsat
206 WHERE (q(:, :, :, 1) < 0.) q(:, :, :, 1) = 1E-10
207 DO l = 1, llm
208 q(:, 1, l, 1) = SUM(aire_2d(:, 1) * q(:, 1, l, 1)) / apoln
209 q(:, jjm + 1, l, 1) &
210 = SUM(aire_2d(:, jjm + 1) * q(:, jjm + 1, l, 1)) / apols
211 ENDDO
212
213 q(:, :, :, 2:4) = 0. ! liquid water, radon and lead
214
215 if (nqmx >= 5) then
216 ! Ozone:
217 call regr_lat_time_coefoz
218 call regr_pr_o3(p3d, q(:, :, :, 5))
219 ! Convert from mole fraction to mass fraction:
220 q(:, :, :, 5) = q(:, :, :, 5) * 48. / 29.
221 end if
222
223 null_array = 0.
224 rugmer = 0.001
225 zmea = pack(zmea_2d, dyn_phy)
226 zstd = pack(zstd_2d, dyn_phy)
227 zsig = pack(zsig_2d, dyn_phy)
228 zgam = pack(zgam_2d, dyn_phy)
229 zthe = pack(zthe_2d, dyn_phy)
230 zpic = pack(zpic_2d, dyn_phy)
231 zval = pack(zval_2d, dyn_phy)
232
233 ! On initialise les sous-surfaces.
234 ! Lecture du fichier glace de terre pour fixer la fraction de terre
235 ! et de glace de terre :
236
237 call nf95_open("landiceref.nc", nf90_nowrite, ncid)
238
239 call nf95_inq_varid(ncid, 'longitude', varid)
240 call nf95_gw_var(ncid, varid, dlon_lic)
241 iml_lic = size(dlon_lic)
242
243 call nf95_inq_varid(ncid, 'latitude', varid)
244 call nf95_gw_var(ncid, varid, dlat_lic)
245 jml_lic = size(dlat_lic)
246
247 call nf95_inq_varid(ncid, 'landice', varid)
248 ALLOCATE(fraclic(iml_lic, jml_lic))
249 call nf95_get_var(ncid, varid, fraclic)
250
251 call nf95_close(ncid)
252
253 ! Interpolation sur la grille T du mod\`ele :
254 PRINT *, 'Dimensions de "landiceref.nc"'
255 print *, "iml_lic = ", iml_lic
256 print *, "jml_lic = ", jml_lic
257
258 ! Si les coordonn\'ees sont en degr\'es, on les transforme :
259 IF (MAXVAL(dlon_lic) > pi) THEN
260 dlon_lic = dlon_lic * pi / 180.
261 ENDIF
262 IF (maxval(dlat_lic) > pi) THEN
263 dlat_lic = dlat_lic * pi/ 180.
264 ENDIF
265
266 flic_tmp(:iim, :) = grille_m(dlon_lic, dlat_lic, fraclic, rlonv(:iim), &
267 rlatu)
268 flic_tmp(iim + 1, :) = flic_tmp(1, :)
269
270 deallocate(dlon_lic, dlat_lic) ! pointers
271
272 ! Passage sur la grille physique
273 pctsrf = 0.
274 pctsrf(:, is_lic) = pack(flic_tmp, dyn_phy)
275 ! Ad\'equation avec le maque terre/mer
276 WHERE (pctsrf(:, is_lic) < EPSFRA) pctsrf(:, is_lic) = 0.
277 WHERE (zmasq < EPSFRA) pctsrf(:, is_lic) = 0.
278 pctsrf(:, is_ter) = zmasq
279 where (zmasq > EPSFRA)
280 where (pctsrf(:, is_lic) >= zmasq)
281 pctsrf(:, is_lic) = zmasq
282 pctsrf(:, is_ter) = 0.
283 elsewhere
284 pctsrf(:, is_ter) = zmasq - pctsrf(:, is_lic)
285 where (pctsrf(:, is_ter) < EPSFRA)
286 pctsrf(:, is_ter) = 0.
287 pctsrf(:, is_lic) = zmasq
288 end where
289 end where
290 end where
291
292 ! Sous-surface oc\'ean et glace de mer (pour d\'emarrer on met glace
293 ! de mer \`a 0) :
294 pctsrf(:, is_oce) = 1. - zmasq
295 WHERE (pctsrf(:, is_oce) < EPSFRA) pctsrf(:, is_oce) = 0.
296
297 ! V\'erification que somme des sous-surfaces vaut 1 :
298 ji = count(abs(sum(pctsrf, dim = 2) - 1.) > EPSFRA)
299 IF (ji /= 0) then
300 PRINT *, 'Bad surface percentages for ', ji, 'points'
301 end IF
302
303 ! Calcul interm\'ediaire :
304 CALL massdair(p3d, masse)
305
306 forall (l = 1:llm)
307 masse(:, 1, l) = SUM(aire_2d(:iim, 1) * masse(:iim, 1, l)) / apoln
308 masse(:, jjm + 1, l) = &
309 SUM(aire_2d(:iim, jjm + 1) * masse(:iim, jjm + 1, l)) / apols
310 END forall
311
312 call iniadvtrac
313 CALL geopot(teta, pk , pks, phis, phi)
314 CALL caldyn0(ucov, vcov, teta, ps, masse, pk, phis, phi, w, pbaru, pbarv)
315 CALL dynredem0(day_ref, phis)
316 CALL dynredem1(vcov, ucov, teta, q, masse, ps, itau = 0)
317
318 ! Initialisations :
319 snsrf = 0.
320 albe(:, is_ter) = 0.08
321 albe(:, is_lic) = 0.6
322 albe(:, is_oce) = 0.5
323 albe(:, is_sic) = 0.6
324 evap = 0.
325 qsolsrf = 150.
326 tsoil = spread(spread(pack(tsol_2d, dyn_phy), 2, nsoilmx), 3, nbsrf)
327 solsw = 165.
328 sollw = -53.
329 t_ancien = 273.15
330 q_ancien = 0.
331 agesno = 0.
332
333 frugs(:, is_oce) = rugmer
334 frugs(:, is_ter) = MAX(1e-5, zstd * zsig / 2)
335 frugs(:, is_lic) = MAX(1e-5, zstd * zsig / 2)
336 frugs(:, is_sic) = 0.001
337 clwcon = 0.
338 rnebcon = 0.
339 ratqs = 0.
340 sig1 = 0.
341 w01 = 0.
342
343 itau_phy = 0
344 nday = 0
345 call phyredem0(lmt_pas = day_step / iphysiq)
346
347 call nf95_inq_varid(ncid_restartphy, "trs", varid)
348 call nf95_put_var(ncid_restartphy, varid, null_array)
349
350 call phyredem(pctsrf, tsoil(:, 1, :), tsoil, tsoil(:, 1, is_oce), &
351 null_array, qsolsrf, pack(qsol_2d, dyn_phy), snsrf, albe, evap, &
352 null_array, null_array, solsw, sollw, null_array, null_array, frugs, &
353 agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, t_ancien, &
354 q_ancien, rnebcon, ratqs, clwcon, null_array, sig1, w01)
355
356 END SUBROUTINE etat0
357
358 end module etat0_mod

  ViewVC Help
Powered by ViewVC 1.1.21