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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.21