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

Annotation of /trunk/dyn3d/etat0.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 13 - (hide annotations)
Fri Jul 25 19:59:34 2008 UTC (15 years, 10 months ago) by guez
Original Path: trunk/libf/dyn3d/etat0.f90
File size: 11935 byte(s)
-- Minor change of behaviour:

"etat0" does not compute "rugsrel" nor "radpas". Deleted arguments
"radpas" and "rugsrel" of "phyredem". Deleted argument "rugsrel" of
"phyetat0". "startphy.nc" does not contain the variable "RUGSREL". In
"physiq", "rugoro" is set to 0 if not "ok_orodr". The whole program
"etat0_lim" does not use "clesphys2".

-- Minor modification of input/output:

Created subroutine "read_clesphys2". Variables of "clesphys2" are read
in "read_clesphys2" instead of "conf_gcm". "printflag" does not print
variables of "clesphys2".

-- Should not change any result at run time:

References to module "numer_rec" instead of individual modules of
"Numer_rec_Lionel".

Deleted argument "clesphy0" of "calfis", "physiq", "conf_gcm",
"leapfrog", "phyetat0". Deleted variable "clesphy0" in
"gcm". "phyetat0" does not modify variables of "clesphys2".

The program unit "gcm" does not modify "itau_phy".

Added some "intent" attributes.

"regr11_lint" does not call "polint".

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

  ViewVC Help
Powered by ViewVC 1.1.21