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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3 - (hide annotations)
Wed Feb 27 13:16:39 2008 UTC (16 years, 3 months ago) by guez
Original Path: trunk/libf/dyn3d/etat0.f90
File size: 12980 byte(s)
Initial import
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 clesphys, only: ok_orodr, nbapp_rad
40     use grid_atob, only: grille_m
41     use grid_change, only: init_dyn_phy, dyn_phy
42     use q_sat_m, only: q_sat
43     use exner_hyb_m, only: exner_hyb
44     use regr_coefoz_m, only: regr_coefoz
45     use advtrac_m, only: iniadvtrac
46     use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, nf90_nowrite, &
47     nf90_get_var, handle_err
48     use pressure_m, only: pls, p3d
49    
50     ! Variables local to the procedure:
51    
52     REAL latfi(klon), lonfi(klon)
53     ! (latitude and longitude of a point of the scalar grid identified
54     ! by a simple index, in °)
55    
56     REAL, dimension(iim + 1, jjm + 1, llm):: uvent, t3d, tpot
57     REAL vvent(iim + 1, jjm, llm)
58    
59     REAL q3d(iim + 1, jjm + 1, llm, nqmx)
60     ! (mass fractions of trace species
61     ! "q3d(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)"
62     ! and pressure level "pls(i, j, l)".)
63    
64     real qsat(iim + 1, jjm + 1, llm) ! mass fraction of saturating water vapor
65     REAL tsol(klon), qsol(klon), sn(klon)
66     REAL tsolsrf(klon, nbsrf), qsolsrf(klon, nbsrf), snsrf(klon, nbsrf)
67     REAL albe(klon, nbsrf), evap(klon, nbsrf)
68     REAL alblw(klon, nbsrf)
69     REAL tsoil(klon, nsoilmx, nbsrf)
70     REAL radsol(klon), rain_fall(klon), snow_fall(klon)
71     REAL solsw(klon), sollw(klon), fder(klon)
72     !IM "slab" ocean
73     REAL tslab(klon)
74     real seaice(klon) ! kg m-2
75     REAL frugs(klon, nbsrf), agesno(klon, nbsrf)
76     REAL rugmer(klon)
77     real, dimension(iim + 1, jjm + 1):: relief, zstd_2d, zsig_2d, zgam_2d
78     real, dimension(iim + 1, jjm + 1):: zthe_2d, zpic_2d, zval_2d
79     real, dimension(iim + 1, jjm + 1):: tsol_2d, psol
80     REAL zmea(klon), zstd(klon)
81     REAL zsig(klon), zgam(klon)
82     REAL zthe(klon)
83     REAL zpic(klon), zval(klon)
84     REAL rugsrel(klon)
85     REAL t_ancien(klon, llm), q_ancien(klon, llm) !
86     REAL run_off_lic_0(klon)
87     real clwcon(klon, llm), rnebcon(klon, llm), ratqs(klon, llm)
88     ! déclarations pour lecture glace de mer
89     INTEGER iml_lic, jml_lic, llm_tmp, ttm_tmp
90     INTEGER itaul(1), fid
91     REAL lev(1), date
92     REAL, ALLOCATABLE:: lon_lic(:, :), lat_lic(:, :)
93     REAL, ALLOCATABLE:: dlon_lic(:), dlat_lic(:)
94     REAL, ALLOCATABLE:: fraclic(:, :) ! fraction land ice
95     REAL flic_tmp(iim + 1, jjm + 1) !fraction land ice temporary
96    
97     INTEGER l, ji
98     INTEGER nq
99    
100     REAL pk(iim + 1, jjm + 1, llm) ! fonction d'Exner aux milieux des couches
101     real pks(iim + 1, jjm + 1)
102    
103     REAL masse(iim + 1, jjm + 1, llm)
104     REAL phi(ip1jmp1, llm)
105     REAL pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
106     REAL w(ip1jmp1, llm)
107     REAL phystep
108     INTEGER radpas
109     integer ncid, varid, ncerr, month
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     zval_2d) ! also compute "masque" and "phis"
137     call init_dyn_phy ! define the mask "dyn_phy" for distinct grid points
138     zmasq = pack(masque, dyn_phy)
139     PRINT *, 'Masque construit'
140    
141     CALL start_init_dyn(tsol_2d, psol) ! also compute "qsol_2d"
142    
143     ! Compute pressure on intermediate levels:
144     forall(l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * psol(:, :)
145     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     uvent(:, :, :) = start_inter_3d('U', rlonv, rlatv, pls)
153     forall (l = 1: llm) uvent(:iim, :, l) = uvent(:iim, :, l) * cu_2d(:iim, :)
154     uvent(iim+1, :, :) = uvent(1, :, :)
155    
156     vvent(:, :, :) = start_inter_3d('V', rlonu, rlatu(:jjm), pls(:, :jjm, :))
157     forall (l = 1: llm) vvent(:iim, :, l) = vvent(:iim, :, l) * cv_2d(:iim, :)
158     vvent(iim + 1, :, :) = vvent(1, :, :)
159    
160     t3d(:, :, :) = start_inter_3d('TEMP', rlonu, rlatv, pls)
161     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     q3d(:, :, :, 1) = 0.01 * start_inter_3d('R', rlonu, rlatv, pls) * qsat
180     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     ! Ozone:
190    
191     ! Compute ozone parameters on the LMDZ grid:
192     call regr_coefoz
193    
194     ! Find the month containing the day number "dayref":
195     month = (dayref - 1) / 30 + 1
196     print *, "month = ", month
197    
198     call nf95_open("coefoz_LMDZ.nc", nf90_nowrite, ncid)
199    
200     ! Get data at the right month from the input file:
201     call nf95_inq_varid(ncid, "r_Mob", varid)
202     ncerr = nf90_get_var(ncid, varid, q3d(:, :, :, 5), &
203     start=(/1, 1, 1, month/))
204     call handle_err("nf90_get_var r_Mob", ncerr)
205    
206     call nf95_close(ncid)
207     ! Latitudes are in increasing order in the input file while
208     ! "rlatu" is in decreasing order so we need to invert order. Also, we
209     ! compute mass fraction from mole fraction:
210     q3d(:, :, :, 5) = q3d(:, jjm+1:1:-1, :, 5) * 48. / 29.
211    
212     tsol(:) = pack(tsol_2d, dyn_phy)
213     qsol(:) = pack(qsol_2d, dyn_phy)
214     sn(:) = 0. ! snow
215     radsol(:) = 0.
216     tslab(:) = 0. ! IM "slab" ocean
217     seaice(:) = 0.
218     rugmer(:) = 0.001
219     zmea(:) = pack(relief, dyn_phy)
220     zstd(:) = pack(zstd_2d, dyn_phy)
221     zsig(:) = pack(zsig_2d, dyn_phy)
222     zgam(:) = pack(zgam_2d, dyn_phy)
223     zthe(:) = pack(zthe_2d, dyn_phy)
224     zpic(:) = pack(zpic_2d, dyn_phy)
225     zval(:) = pack(zval_2d, dyn_phy)
226    
227     rugsrel(:) = 0.
228     IF (ok_orodr) rugsrel(:) = MAX(1.e-05, zstd(:) * zsig(:) / 2)
229    
230     ! On initialise les sous-surfaces:
231     ! Lecture du fichier glace de terre pour fixer la fraction de terre
232     ! et de glace de terre:
233     CALL flininfo("landiceref.nc", iml_lic, jml_lic, llm_tmp, &
234     ttm_tmp, fid)
235     ALLOCATE(lat_lic(iml_lic, jml_lic))
236     ALLOCATE(lon_lic(iml_lic, jml_lic))
237     ALLOCATE(dlon_lic(iml_lic))
238     ALLOCATE(dlat_lic(jml_lic))
239     ALLOCATE(fraclic(iml_lic, jml_lic))
240     CALL flinopen_nozoom("landiceref.nc", iml_lic, jml_lic, &
241     llm_tmp, lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, &
242     fid)
243     CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp &
244     , 1, 1, fraclic)
245     CALL flinclo(fid)
246    
247     ! Interpolation sur la grille T du modèle :
248     PRINT *, 'Dimensions de "landice"'
249     print *, "iml_lic = ", iml_lic
250     print *, "jml_lic = ", jml_lic
251    
252     ! Si les coordonnées sont en degrés, on les transforme :
253     IF (MAXVAL( lon_lic(:, :) ) > pi) THEN
254     lon_lic(:, :) = lon_lic(:, :) * pi / 180.
255     ENDIF
256     IF (maxval( lat_lic(:, :) ) > pi) THEN
257     lat_lic(:, :) = lat_lic(:, :) * pi/ 180.
258     ENDIF
259    
260     dlon_lic(:) = lon_lic(:, 1)
261     dlat_lic(:) = lat_lic(1, :)
262    
263     flic_tmp(:iim, :) = grille_m(dlon_lic, dlat_lic, fraclic, rlonv(:iim), &
264     rlatu)
265     flic_tmp(iim + 1, :) = flic_tmp(1, :)
266    
267     ! Passage sur la grille physique
268     pctsrf(:, :)=0.
269     pctsrf(:, is_lic) = pack(flic_tmp, dyn_phy)
270     ! Adéquation avec le maque terre/mer
271     WHERE (pctsrf(:, is_lic) < EPSFRA ) pctsrf(:, is_lic) = 0.
272     WHERE (zmasq(:) < EPSFRA) pctsrf(:, is_lic) = 0.
273     pctsrf(:, is_ter) = zmasq(:)
274     where (zmasq(:) > EPSFRA)
275     where (pctsrf(:, is_lic) >= zmasq(:))
276     pctsrf(:, is_lic) = zmasq(:)
277     pctsrf(:, is_ter) = 0.
278     elsewhere
279     pctsrf(:, is_ter) = zmasq(:) - pctsrf(:, is_lic)
280     where (pctsrf(:, is_ter) < EPSFRA)
281     pctsrf(:, is_ter) = 0.
282     pctsrf(:, is_lic) = zmasq(:)
283     end where
284     end where
285     end where
286    
287     ! Sous-surface océan et glace de mer (pour démarrer on met glace
288     ! de mer à 0) :
289     pctsrf(:, is_oce) = 1. - zmasq(:)
290     WHERE (pctsrf(:, is_oce) < EPSFRA) pctsrf(:, is_oce) = 0.
291    
292     ! Vérification que somme des sous-surfaces vaut 1:
293     ji = count(abs(sum(pctsrf(:, :), dim = 2) - 1. ) > EPSFRA)
294     IF (ji /= 0) PRINT *, 'Problème répartition sous maille pour', ji, 'points'
295    
296     ! Calcul intermédiaire:
297     CALL massdair(p3d, masse)
298    
299     print *, 'ALPHAX = ', alphax
300    
301     forall (l = 1:llm)
302     masse(:, 1, l) = SUM(aire_2d(:iim, 1) * masse(:iim, 1, l)) / apoln
303     masse(:, jjm + 1, l) = &
304     SUM(aire_2d(:iim, jjm + 1) * masse(:iim, jjm + 1, l)) / apols
305     END forall
306    
307     ! Initialisation pour traceurs:
308     call iniadvtrac(nq)
309     ! Ecriture:
310     CALL inidissip(lstardis, nitergdiv, nitergrot, niterh, tetagdiv, &
311     tetagrot, tetatemp)
312     itau_dyn = 0
313     itau_phy = 0
314     day_ref = dayref
315     annee_ref = anneeref
316    
317     CALL geopot(ip1jmp1, tpot, pk , pks, phis , phi)
318     CALL caldyn0(0, uvent, vvent, tpot, psol, masse, pk, phis, phi, w, &
319     pbaru, pbarv, 0)
320     CALL dynredem0("start.nc", dayref, phis, nqmx)
321     CALL dynredem1("start.nc", 0., vvent, uvent, tpot, q3d, nqmx, masse, psol)
322    
323     ! Ecriture état initial physique:
324     print *, 'dtvr = ', dtvr
325     print *, "iphysiq = ", iphysiq
326     print *, "nbapp_rad = ", nbapp_rad
327     phystep = dtvr * REAL(iphysiq)
328     radpas = NINT (86400./phystep/ nbapp_rad)
329     print *, 'phystep = ', phystep
330     print *, "radpas = ", radpas
331    
332     ! Initialisations :
333     tsolsrf(:, is_ter) = tsol
334     tsolsrf(:, is_lic) = tsol
335     tsolsrf(:, is_oce) = tsol
336     tsolsrf(:, is_sic) = tsol
337     snsrf(:, is_ter) = sn
338     snsrf(:, is_lic) = sn
339     snsrf(:, is_oce) = sn
340     snsrf(:, is_sic) = sn
341     albe(:, is_ter) = 0.08
342     albe(:, is_lic) = 0.6
343     albe(:, is_oce) = 0.5
344     albe(:, is_sic) = 0.6
345     alblw = albe
346     evap(:, :) = 0.
347     qsolsrf(:, is_ter) = 150.
348     qsolsrf(:, is_lic) = 150.
349     qsolsrf(:, is_oce) = 150.
350     qsolsrf(:, is_sic) = 150.
351     tsoil = spread(spread(tsol, 2, nsoilmx), 3, nbsrf)
352     rain_fall = 0.
353     snow_fall = 0.
354     solsw = 165.
355     sollw = -53.
356     t_ancien = 273.15
357     q_ancien = 0.
358     agesno = 0.
359     !IM "slab" ocean
360     tslab(:) = tsolsrf(:, is_oce)
361     seaice = 0.
362    
363     frugs(:, is_oce) = rugmer(:)
364     frugs(:, is_ter) = MAX(1.e-05, zstd(:) * zsig(:) / 2)
365     frugs(:, is_lic) = MAX(1.e-05, zstd(:) * zsig(:) / 2)
366     frugs(:, is_sic) = 0.001
367     fder = 0.
368     clwcon = 0.
369     rnebcon = 0.
370     ratqs = 0.
371     run_off_lic_0 = 0.
372    
373     call phyredem("startphy.nc", phystep, radpas, latfi, lonfi, pctsrf, &
374     tsolsrf, tsoil, tslab, seaice, qsolsrf, qsol, snsrf, albe, alblw, &
375     evap, rain_fall, snow_fall, solsw, sollw, fder, radsol, frugs, &
376     agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, rugsrel, &
377     t_ancien, q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0)
378     CALL histclo
379    
380     END SUBROUTINE etat0
381    
382     end module etat0_mod

  ViewVC Help
Powered by ViewVC 1.1.21