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

Annotation of /trunk/Sources/dyn3d/etat0.f

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.21