/[lmdze]/trunk/libf/dyn3d/etat0.f90
ViewVC logotype

Annotation of /trunk/libf/dyn3d/etat0.f90

Parent Directory Parent Directory | Revision Log Revision Log


Revision 27 - (hide annotations)
Thu Mar 25 14:29:07 2010 UTC (14 years, 1 month ago) by guez
File size: 11920 byte(s)
"dyn3d" and "filtrez" do not contain any included file so make rules
have been updated.

"comdissip.f90" was useless, removed it.

"dynredem0" wrote undefined value in "controle(31)", that was
overwritten by "dynredem1". Now "dynredem0" just writes 0 to
"controle(31)".

Removed arguments of "inidissip". "inidissip" now accesses the
variables by use association.

In program "etat0_lim", "itaufin" is not defined so "dynredem1" wrote
undefined value to "controle(31)". Added argument "itau" of
"dynredem1" to correct that.

"itaufin" does not need to be a module variable (of "temps"), made it
a local variable of "leapfrog".

Removed calls to "diagedyn" from "leapfrog".

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     ! position "i")
11 guez 3
12     private nbsrf, klon
13    
14     contains
15    
16     SUBROUTINE etat0
17    
18     ! From "etat0_netcdf.F", version 1.3 2005/05/25 13:10:09
19    
20 guez 15 ! This subroutine creates "mask".
21 guez 3
22 guez 27 use caldyn0_m, only: caldyn0
23 guez 3 use comconst, only: dtvr, daysec, cpp, kappa, pi
24 guez 27 use comgeom, only: rlatu, rlonv, rlonu, rlatv, aire_2d, apoln, apols, &
25     cu_2d, cv_2d
26 guez 3 use comvert, only: ap, bp, preff, pa
27 guez 27 use conf_gcm_m, only: day_step, iphysiq, dayref, anneeref
28     use dimens_m, only: iim, jjm, llm, nqmx
29 guez 3 use dimphy, only: zmasq
30     use dimsoil, only: nsoilmx
31 guez 27 use dynredem0_m, only: dynredem0
32     use dynredem1_m, only: dynredem1
33     use exner_hyb_m, only: exner_hyb
34 guez 3 use grid_atob, only: grille_m
35     use grid_change, only: init_dyn_phy, dyn_phy
36 guez 27 use indicesol, only: is_oce, is_sic, is_ter, is_lic, epsfra
37 guez 18 use iniadvtrac_m, only: iniadvtrac
38 guez 27 use inidissip_m, only: inidissip
39     use inigeom_m, only: inigeom
40     USE ioipsl, only: flinget, flinclo, flinopen_nozoom, flininfo, histclo
41     use paramet_m, only: ip1jm, ip1jmp1
42     use phyredem_m, only: phyredem
43 guez 10 use pressure_var, only: pls, p3d
44 guez 27 use q_sat_m, only: q_sat
45 guez 7 use regr_lat_time_coefoz_m, only: regr_lat_time_coefoz
46     use regr_pr_o3_m, only: regr_pr_o3
47 guez 27 use serre, only: alphax
48     USE start_init_orog_m, only: start_init_orog, mask, phis
49     use start_init_phys_m, only: qsol_2d
50     use startdyn, only: start_inter_3d, start_init_dyn
51     use temps, only: itau_phy, annee_ref, day_ref, dt
52 guez 3
53     ! Variables local to the procedure:
54    
55     REAL latfi(klon), lonfi(klon)
56     ! (latitude and longitude of a point of the scalar grid identified
57     ! by a simple index, in °)
58    
59     REAL, dimension(iim + 1, jjm + 1, llm):: uvent, t3d, tpot
60     REAL vvent(iim + 1, jjm, llm)
61    
62     REAL q3d(iim + 1, jjm + 1, llm, nqmx)
63     ! (mass fractions of trace species
64     ! "q3d(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)"
65     ! and pressure level "pls(i, j, l)".)
66    
67     real qsat(iim + 1, jjm + 1, llm) ! mass fraction of saturating water vapor
68     REAL tsol(klon), qsol(klon), sn(klon)
69     REAL tsolsrf(klon, nbsrf), qsolsrf(klon, nbsrf), snsrf(klon, nbsrf)
70     REAL albe(klon, nbsrf), evap(klon, nbsrf)
71     REAL alblw(klon, nbsrf)
72     REAL tsoil(klon, nsoilmx, nbsrf)
73     REAL radsol(klon), rain_fall(klon), snow_fall(klon)
74     REAL solsw(klon), sollw(klon), fder(klon)
75     !IM "slab" ocean
76     REAL tslab(klon)
77     real seaice(klon) ! kg m-2
78     REAL frugs(klon, nbsrf), agesno(klon, nbsrf)
79     REAL rugmer(klon)
80     real, dimension(iim + 1, jjm + 1):: relief, zstd_2d, zsig_2d, zgam_2d
81     real, dimension(iim + 1, jjm + 1):: zthe_2d, zpic_2d, zval_2d
82     real, dimension(iim + 1, jjm + 1):: tsol_2d, psol
83     REAL zmea(klon), zstd(klon)
84     REAL zsig(klon), zgam(klon)
85     REAL zthe(klon)
86     REAL zpic(klon), zval(klon)
87     REAL t_ancien(klon, llm), q_ancien(klon, llm) !
88     REAL run_off_lic_0(klon)
89     real clwcon(klon, llm), rnebcon(klon, llm), ratqs(klon, llm)
90     ! déclarations pour lecture glace de mer
91     INTEGER iml_lic, jml_lic, llm_tmp, ttm_tmp
92     INTEGER itaul(1), fid
93     REAL lev(1), date
94     REAL, ALLOCATABLE:: lon_lic(:, :), lat_lic(:, :)
95     REAL, ALLOCATABLE:: dlon_lic(:), dlat_lic(:)
96     REAL, ALLOCATABLE:: fraclic(:, :) ! fraction land ice
97     REAL flic_tmp(iim + 1, jjm + 1) !fraction land ice temporary
98    
99     INTEGER l, ji
100    
101     REAL pk(iim + 1, jjm + 1, llm) ! fonction d'Exner aux milieux des couches
102     real pks(iim + 1, jjm + 1)
103    
104     REAL masse(iim + 1, jjm + 1, llm)
105     REAL phi(ip1jmp1, llm)
106     REAL pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
107     REAL w(ip1jmp1, llm)
108     REAL phystep
109    
110     !---------------------------------
111    
112     print *, "Call sequence information: etat0"
113    
114     ! Construct a grid:
115    
116     dtvr = daysec / real(day_step)
117     print *, 'dtvr = ', dtvr
118    
119     pa = 5e4
120     CALL iniconst
121     CALL inigeom
122     CALL inifilr
123    
124     latfi(1) = 90.
125     latfi(2:klon-1) = pack(spread(rlatu(2:jjm), 1, iim), .true.) * 180. / pi
126     ! (with conversion to degrees)
127     latfi(klon) = - 90.
128    
129     lonfi(1) = 0.
130     lonfi(2:klon-1) = pack(spread(rlonv(:iim), 2, jjm - 1), .true.) * 180. / pi
131     ! (with conversion to degrees)
132     lonfi(klon) = 0.
133    
134     call start_init_orog(relief, zstd_2d, zsig_2d, zgam_2d, zthe_2d, zpic_2d, &
135 guez 15 zval_2d) ! also compute "mask" and "phis"
136 guez 3 call init_dyn_phy ! define the mask "dyn_phy" for distinct grid points
137 guez 15 zmasq = pack(mask, dyn_phy)
138 guez 3 PRINT *, 'Masque construit'
139    
140     CALL start_init_dyn(tsol_2d, psol) ! also compute "qsol_2d"
141    
142     ! Compute pressure on intermediate levels:
143 guez 15 forall(l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * psol
144 guez 3 CALL exner_hyb(psol, p3d, pks, pk)
145     IF (MINVAL(pk) == MAXVAL(pk)) stop '"pk" should not be constant'
146    
147     pls(:, :, :) = preff * (pk(:, :, :) / cpp)**(1. / kappa)
148     PRINT *, "minval(pls(:, :, :)) = ", minval(pls(:, :, :))
149     print *, "maxval(pls(:, :, :)) = ", maxval(pls(:, :, :))
150    
151 guez 23 call start_inter_3d('U', rlonv, rlatv, pls, uvent)
152 guez 3 forall (l = 1: llm) uvent(:iim, :, l) = uvent(:iim, :, l) * cu_2d(:iim, :)
153     uvent(iim+1, :, :) = uvent(1, :, :)
154    
155 guez 23 call start_inter_3d('V', rlonu, rlatu(:jjm), pls(:, :jjm, :), vvent)
156 guez 3 forall (l = 1: llm) vvent(:iim, :, l) = vvent(:iim, :, l) * cv_2d(:iim, :)
157     vvent(iim + 1, :, :) = vvent(1, :, :)
158    
159 guez 23 call start_inter_3d('TEMP', rlonu, rlatv, pls, t3d)
160 guez 3 PRINT *, 'minval(t3d(:, :, :)) = ', minval(t3d(:, :, :))
161     print *, "maxval(t3d(:, :, :)) = ", maxval(t3d(:, :, :))
162    
163     tpot(:iim, :, :) = t3d(:iim, :, :) * cpp / pk(:iim, :, :)
164     tpot(iim + 1, :, :) = tpot(1, :, :)
165     DO l=1, llm
166     tpot(:, 1, l) = SUM(aire_2d(:, 1) * tpot(:, 1, l)) / apoln
167     tpot(:, jjm + 1, l) = SUM(aire_2d(:, jjm + 1) * tpot(:, jjm + 1, l)) &
168     / apols
169     ENDDO
170    
171     ! Calcul de l'humidité à saturation :
172     qsat(:, :, :) = q_sat(t3d, pls)
173     PRINT *, "minval(qsat(:, :, :)) = ", minval(qsat(:, :, :))
174     print *, "maxval(qsat(:, :, :)) = ", maxval(qsat(:, :, :))
175     IF (MINVAL(qsat) == MAXVAL(qsat)) stop '"qsat" should not be constant'
176    
177     ! Water vapor:
178 guez 23 call start_inter_3d('R', rlonu, rlatv, pls, q3d(:, :, :, 1))
179     q3d(:, :, :, 1) = 0.01 * q3d(:, :, :, 1) * qsat
180 guez 3 WHERE (q3d(:, :, :, 1) < 0.) q3d(:, :, :, 1) = 1E-10
181     DO l = 1, llm
182     q3d(:, 1, l, 1) = SUM(aire_2d(:, 1) * q3d(:, 1, l, 1)) / apoln
183     q3d(:, jjm + 1, l, 1) &
184     = SUM(aire_2d(:, jjm + 1) * q3d(:, jjm + 1, l, 1)) / apols
185     ENDDO
186    
187     q3d(:, :, :, 2:4) = 0. ! liquid water, radon and lead
188    
189 guez 5 if (nqmx >= 5) then
190     ! Ozone:
191 guez 7 call regr_lat_time_coefoz
192     call regr_pr_o3(q3d(:, :, :, 5))
193     ! Convert from mole fraction to mass fraction:
194     q3d(:, :, :, 5) = q3d(:, :, :, 5) * 48. / 29.
195 guez 5 end if
196 guez 3
197 guez 13 tsol = pack(tsol_2d, dyn_phy)
198     qsol = pack(qsol_2d, dyn_phy)
199     sn = 0. ! snow
200     radsol = 0.
201     tslab = 0. ! IM "slab" ocean
202     seaice = 0.
203     rugmer = 0.001
204     zmea = pack(relief, dyn_phy)
205     zstd = pack(zstd_2d, dyn_phy)
206     zsig = pack(zsig_2d, dyn_phy)
207     zgam = pack(zgam_2d, dyn_phy)
208     zthe = pack(zthe_2d, dyn_phy)
209     zpic = pack(zpic_2d, dyn_phy)
210     zval = pack(zval_2d, dyn_phy)
211 guez 3
212     ! On initialise les sous-surfaces:
213     ! Lecture du fichier glace de terre pour fixer la fraction de terre
214     ! et de glace de terre:
215     CALL flininfo("landiceref.nc", iml_lic, jml_lic, llm_tmp, &
216     ttm_tmp, fid)
217     ALLOCATE(lat_lic(iml_lic, jml_lic))
218     ALLOCATE(lon_lic(iml_lic, jml_lic))
219     ALLOCATE(dlon_lic(iml_lic))
220     ALLOCATE(dlat_lic(jml_lic))
221     ALLOCATE(fraclic(iml_lic, jml_lic))
222     CALL flinopen_nozoom("landiceref.nc", iml_lic, jml_lic, &
223     llm_tmp, lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, &
224     fid)
225     CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp &
226     , 1, 1, fraclic)
227     CALL flinclo(fid)
228    
229     ! Interpolation sur la grille T du modèle :
230     PRINT *, 'Dimensions de "landice"'
231     print *, "iml_lic = ", iml_lic
232     print *, "jml_lic = ", jml_lic
233    
234     ! Si les coordonnées sont en degrés, on les transforme :
235 guez 15 IF (MAXVAL( lon_lic ) > pi) THEN
236     lon_lic = lon_lic * pi / 180.
237 guez 3 ENDIF
238 guez 15 IF (maxval( lat_lic ) > pi) THEN
239     lat_lic = lat_lic * pi/ 180.
240 guez 3 ENDIF
241    
242 guez 13 dlon_lic = lon_lic(:, 1)
243     dlat_lic = lat_lic(1, :)
244 guez 3
245     flic_tmp(:iim, :) = grille_m(dlon_lic, dlat_lic, fraclic, rlonv(:iim), &
246     rlatu)
247     flic_tmp(iim + 1, :) = flic_tmp(1, :)
248    
249     ! Passage sur la grille physique
250 guez 15 pctsrf = 0.
251 guez 3 pctsrf(:, is_lic) = pack(flic_tmp, dyn_phy)
252     ! Adéquation avec le maque terre/mer
253     WHERE (pctsrf(:, is_lic) < EPSFRA ) pctsrf(:, is_lic) = 0.
254 guez 13 WHERE (zmasq < EPSFRA) pctsrf(:, is_lic) = 0.
255     pctsrf(:, is_ter) = zmasq
256     where (zmasq > EPSFRA)
257     where (pctsrf(:, is_lic) >= zmasq)
258     pctsrf(:, is_lic) = zmasq
259 guez 3 pctsrf(:, is_ter) = 0.
260     elsewhere
261 guez 13 pctsrf(:, is_ter) = zmasq - pctsrf(:, is_lic)
262 guez 3 where (pctsrf(:, is_ter) < EPSFRA)
263     pctsrf(:, is_ter) = 0.
264 guez 13 pctsrf(:, is_lic) = zmasq
265 guez 3 end where
266     end where
267     end where
268    
269     ! Sous-surface océan et glace de mer (pour démarrer on met glace
270     ! de mer à 0) :
271 guez 13 pctsrf(:, is_oce) = 1. - zmasq
272 guez 3 WHERE (pctsrf(:, is_oce) < EPSFRA) pctsrf(:, is_oce) = 0.
273    
274     ! Vérification que somme des sous-surfaces vaut 1:
275 guez 15 ji = count(abs(sum(pctsrf, dim = 2) - 1.) > EPSFRA)
276     IF (ji /= 0) then
277     PRINT *, 'Problème répartition sous maille pour ', ji, 'points'
278     end IF
279 guez 3
280     ! Calcul intermédiaire:
281     CALL massdair(p3d, masse)
282    
283     print *, 'ALPHAX = ', alphax
284    
285     forall (l = 1:llm)
286     masse(:, 1, l) = SUM(aire_2d(:iim, 1) * masse(:iim, 1, l)) / apoln
287     masse(:, jjm + 1, l) = &
288     SUM(aire_2d(:iim, jjm + 1) * masse(:iim, jjm + 1, l)) / apols
289     END forall
290    
291     ! Initialisation pour traceurs:
292 guez 5 call iniadvtrac
293 guez 27 CALL inidissip
294 guez 3 itau_phy = 0
295     day_ref = dayref
296     annee_ref = anneeref
297    
298 guez 22 CALL geopot(ip1jmp1, tpot, pk , pks, phis, phi)
299 guez 23 CALL caldyn0(uvent, vvent, tpot, psol, masse, pk, phis, phi, w, pbaru, &
300     pbarv)
301 guez 5 CALL dynredem0("start.nc", dayref, phis)
302 guez 27 CALL dynredem1("start.nc", vvent, uvent, tpot, q3d, masse, psol, itau=0)
303 guez 3
304     ! Ecriture état initial physique:
305     print *, "iphysiq = ", iphysiq
306     phystep = dtvr * REAL(iphysiq)
307     print *, 'phystep = ', phystep
308    
309     ! Initialisations :
310     tsolsrf(:, is_ter) = tsol
311     tsolsrf(:, is_lic) = tsol
312     tsolsrf(:, is_oce) = tsol
313     tsolsrf(:, is_sic) = tsol
314     snsrf(:, is_ter) = sn
315     snsrf(:, is_lic) = sn
316     snsrf(:, is_oce) = sn
317     snsrf(:, is_sic) = sn
318     albe(:, is_ter) = 0.08
319     albe(:, is_lic) = 0.6
320     albe(:, is_oce) = 0.5
321     albe(:, is_sic) = 0.6
322     alblw = albe
323 guez 15 evap = 0.
324 guez 3 qsolsrf(:, is_ter) = 150.
325     qsolsrf(:, is_lic) = 150.
326     qsolsrf(:, is_oce) = 150.
327     qsolsrf(:, is_sic) = 150.
328     tsoil = spread(spread(tsol, 2, nsoilmx), 3, nbsrf)
329     rain_fall = 0.
330     snow_fall = 0.
331     solsw = 165.
332     sollw = -53.
333     t_ancien = 273.15
334     q_ancien = 0.
335     agesno = 0.
336     !IM "slab" ocean
337 guez 13 tslab = tsolsrf(:, is_oce)
338 guez 3 seaice = 0.
339    
340 guez 13 frugs(:, is_oce) = rugmer
341     frugs(:, is_ter) = MAX(1.e-05, zstd * zsig / 2)
342     frugs(:, is_lic) = MAX(1.e-05, zstd * zsig / 2)
343 guez 3 frugs(:, is_sic) = 0.001
344     fder = 0.
345     clwcon = 0.
346     rnebcon = 0.
347     ratqs = 0.
348     run_off_lic_0 = 0.
349    
350 guez 13 call phyredem("startphy.nc", latfi, lonfi, pctsrf, &
351 guez 3 tsolsrf, tsoil, tslab, seaice, qsolsrf, qsol, snsrf, albe, alblw, &
352     evap, rain_fall, snow_fall, solsw, sollw, fder, radsol, frugs, &
353 guez 13 agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, &
354 guez 3 t_ancien, q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0)
355     CALL histclo
356    
357     END SUBROUTINE etat0
358    
359     end module etat0_mod

  ViewVC Help
Powered by ViewVC 1.1.21