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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 25 - (hide annotations)
Fri Mar 5 16:43:45 2010 UTC (14 years, 2 months ago) by guez
Original Path: trunk/libf/dyn3d/etat0.f90
File size: 12069 byte(s)
Simplified "etat0_lim.sh" and "gcm.sh" because the full versions
depended on personal arrangements for directories and machines.

Translated included files into modules. Encapsulated procedures into modules.

Moved variables from module "comgeom" to local variables of
"inigeom". Deleted some unused variables in "comgeom".

Moved variable "day_ini" from module "temps" to module "dynetat0_m".

Removed useless test on variable "time" and useless "close" statement
in procedure "leapfrog".

Removed useless call to "inigeom" in procedure "limit".

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 3
54     ! Variables local to the procedure:
55    
56     REAL latfi(klon), lonfi(klon)
57     ! (latitude and longitude of a point of the scalar grid identified
58     ! by a simple index, in °)
59    
60     REAL, dimension(iim + 1, jjm + 1, llm):: uvent, t3d, tpot
61     REAL vvent(iim + 1, jjm, llm)
62    
63     REAL q3d(iim + 1, jjm + 1, llm, nqmx)
64     ! (mass fractions of trace species
65     ! "q3d(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)"
66     ! and pressure level "pls(i, j, l)".)
67    
68     real qsat(iim + 1, jjm + 1, llm) ! mass fraction of saturating water vapor
69     REAL tsol(klon), qsol(klon), sn(klon)
70     REAL tsolsrf(klon, nbsrf), qsolsrf(klon, nbsrf), snsrf(klon, nbsrf)
71     REAL albe(klon, nbsrf), evap(klon, nbsrf)
72     REAL alblw(klon, nbsrf)
73     REAL tsoil(klon, nsoilmx, nbsrf)
74     REAL radsol(klon), rain_fall(klon), snow_fall(klon)
75     REAL solsw(klon), sollw(klon), fder(klon)
76     !IM "slab" ocean
77     REAL tslab(klon)
78     real seaice(klon) ! kg m-2
79     REAL frugs(klon, nbsrf), agesno(klon, nbsrf)
80     REAL rugmer(klon)
81     real, dimension(iim + 1, jjm + 1):: relief, zstd_2d, zsig_2d, zgam_2d
82     real, dimension(iim + 1, jjm + 1):: zthe_2d, zpic_2d, zval_2d
83     real, dimension(iim + 1, jjm + 1):: tsol_2d, psol
84     REAL zmea(klon), zstd(klon)
85     REAL zsig(klon), zgam(klon)
86     REAL zthe(klon)
87     REAL zpic(klon), zval(klon)
88     REAL t_ancien(klon, llm), q_ancien(klon, llm) !
89     REAL run_off_lic_0(klon)
90     real clwcon(klon, llm), rnebcon(klon, llm), ratqs(klon, llm)
91     ! déclarations pour lecture glace de mer
92     INTEGER iml_lic, jml_lic, llm_tmp, ttm_tmp
93     INTEGER itaul(1), fid
94     REAL lev(1), date
95     REAL, ALLOCATABLE:: lon_lic(:, :), lat_lic(:, :)
96     REAL, ALLOCATABLE:: dlon_lic(:), dlat_lic(:)
97     REAL, ALLOCATABLE:: fraclic(:, :) ! fraction land ice
98     REAL flic_tmp(iim + 1, jjm + 1) !fraction land ice temporary
99    
100     INTEGER l, ji
101    
102     REAL pk(iim + 1, jjm + 1, llm) ! fonction d'Exner aux milieux des couches
103     real pks(iim + 1, jjm + 1)
104    
105     REAL masse(iim + 1, jjm + 1, llm)
106     REAL phi(ip1jmp1, llm)
107     REAL pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
108     REAL w(ip1jmp1, llm)
109     REAL phystep
110    
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     llm_tmp, lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, &
225     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 3 ! Ecriture:
295     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