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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.21