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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 43 - (hide annotations)
Fri Apr 8 12:43:31 2011 UTC (13 years, 1 month ago) by guez
Original Path: trunk/libf/dyn3d/etat0.f90
File size: 11917 byte(s)
"start_init_phys" is now called directly by "etat0" instead of through
"start_init_dyn". "qsol_2d" is no longer a variable of module
"start_init_phys_m", it is an argument of
"start_init_phys". "start_init_dyn" now receives "tsol_2d" from
"etat0".

Split file "vlspltqs.f" into "vlspltqs.f90", "vlxqs.f90" and
""vlyqs.f90".

In "start_init_orog", replaced calls to "flin*" by calls to NetCDF95.

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

  ViewVC Help
Powered by ViewVC 1.1.21