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

Annotation of /trunk/dyn3d/etat0.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 39 - (hide annotations)
Tue Jan 25 15:11:05 2011 UTC (13 years, 4 months ago) by guez
Original Path: trunk/libf/dyn3d/etat0.f90
File size: 11851 byte(s)
"pi" comes from "nr_util". Removed subroutine "initialize" in module
"comconst".

Copied the content of "fxy_sin.h" into "fxysinus", instead of getting
it from an "include" line. Removed file "fxy_sin.h".

"ps" has rank 2 in "gcm" and "dynetat0".

Assumed-shape for argument "q" of "integrd".

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

  ViewVC Help
Powered by ViewVC 1.1.21