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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 28 - (hide annotations)
Fri Mar 26 18:33:04 2010 UTC (14 years, 1 month ago) by guez
File size: 11934 byte(s)
Removed unused "diagedyn.f" and "undefSTD.f".

In "etat0", the variable "dt" of module "temps" was defined from
"landicered.nc", which was meaningless and useless. Replaced "dt" by a
local trash variable.

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

  ViewVC Help
Powered by ViewVC 1.1.21