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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 23 - (hide annotations)
Mon Dec 14 15:25:16 2009 UTC (14 years, 5 months ago) by guez
Original Path: trunk/libf/dyn3d/etat0.f90
File size: 12065 byte(s)
Split "orografi.f": one file for each procedure. Put the created files
in new directory "Orography".

Removed argument "vcov" of procedure "sortvarc". Removed arguments
"itau" and "time" of procedure "caldyn0". Removed arguments "itau",
"time" and "vcov" of procedure "sortvarc0".

Removed argument "time" of procedure "dynredem1". Removed NetCDF
variable "temps" in files "start.nc" and "restart.nc", because its
value is always 0.

Removed argument "nq" of procedures "iniadvtrac" and "leapfrog". The
number of "tracers read in "traceur.def" must now be equal to "nqmx",
or "nqmx" must equal 4 if there is no file "traceur.def". Replaced
variable "nq" by constant "nqmx" in "leapfrog".

NetCDF variable for ozone field in "coefoz.nc" must now be called
"tro3" instead of "r".

Fixed bug in "zenang".

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