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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 191 - (hide annotations)
Mon May 9 19:56:28 2016 UTC (8 years ago) by guez
File size: 11614 byte(s)
Extracted the call to read_comdissnew out of conf_gcm.

Made ok_instan a variable of module clesphys, itau_phy a variable of
module phyetat0_m, nid_ins a variable of module ini_histins_m, itap a
variable of new module time_phylmdz, so that histwrite_phy can be
called from any procedure without the need to cascade those variables
into that procedure. Made itau_w a variable of module time_phylmdz so
that it is computed only once per time step of physics.

Extracted variables of module clesphys which were in namelist
conf_phys_nml into their own namelist, clesphys_nml, and created
procedure read_clesphys reading clesphys_nml, to avoid side effect.

No need for double precision in procedure getso4fromfile. Assume there
is a single variable for the whole year in the NetCDF file instead of
one variable per month.

Created generic procedure histwrite_phy and removed procedure
write_histins, following LMDZ. histwrite_phy has only two arguments,
can be called from anywhere, and should manage the logic of writing or
not writing into various history files with various operations. So the
test on ok_instan goes inside histwrite_phy.

Test for raz_date in phyetat0 instead of physiq to avoid side effect.

Created procedure increment_itap to avoid side effect.

Removed unnecessary differences between procedures readsulfate and
readsulfate_pi.

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

  ViewVC Help
Powered by ViewVC 1.1.21