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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 26 - (hide annotations)
Tue Mar 9 15:27:15 2010 UTC (14 years, 2 months ago) by guez
File size: 12090 byte(s)
Moved variable "dtdiss" from module "comconst", variable "idissip"
from module "conf_gcm_m" and all variables from module "comdissipn" to
module "inidissip_m". "inidissip" creates file
"inidissip.csv". "idissip" is no longer read from a namelist. Removed
useless computation of "dtdiss" in procedure "iniconst".

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