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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 7 - (hide annotations)
Mon Mar 31 12:24:17 2008 UTC (16 years, 1 month ago) by guez
File size: 12343 byte(s)
This revision is not in working order. Pending some moving of files.

Important changes. In the program "etat0_lim": ozone coefficients from
Mobidic are regridded in time instead of pressure ; consequences in
"etat0". In the program "gcm", ozone coefficients from Mobidic are
read once per day only for the current day and regridded in pressure ;
consequences in "o3_chem_m", "regr_pr_coefoz", "phytrac" and
"regr_pr_comb_coefoz_m".

NetCDF95 is a library and does not export NetCDF.

New variables "nag_gl_options", "nag_fcalls_options" and
"nag_cross_options" in "nag_tools.mk".

"check_coefoz.jnl" rewritten entirely for new version of
"coefoz_LMDZ.nc".

Target "obj_etat0_lim" moved from "GNUmakefile" to "nag_rules.mk".

Added some "intent" attributes in "calfis", "clmain", "clqh",
"cltrac", "cltracrn", "cvltr", "ini_undefSTD", "moy_undefSTD",
"nflxtr", "phystokenc", "phytrac", "readsulfate", "readsulfate_preind"
and "undefSTD".

In "dynetat0", "dynredem0" and "gcm", "phis" has rank 2 instead of
1. "phis" has assumed shape in "dynredem0".

Added module containing "dynredem0". Changed some calls with NetCDF
Fortran 77 interface to calls with NetCDF95 interface.

Replaced calls to "ssum" by calls to "sum" in "inigeom".

In "make.sh", new option "-c" to change compiler.

In "aaam_bud", argument "rjour" deleted.

In "physiq": renamed some variables; deleted variable "xjour".

In "phytrac": renamed some variables; new argument "lmt_pas".

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 advtrac_m, only: iniadvtrac
45     use pressure_m, only: pls, p3d
46 guez 7 use dynredem0_m, only: dynredem0
47     use regr_lat_time_coefoz_m, only: regr_lat_time_coefoz
48     use regr_pr_o3_m, only: regr_pr_o3
49 guez 3
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    
99     REAL pk(iim + 1, jjm + 1, llm) ! fonction d'Exner aux milieux des couches
100     real pks(iim + 1, jjm + 1)
101    
102     REAL masse(iim + 1, jjm + 1, llm)
103     REAL phi(ip1jmp1, llm)
104     REAL pbaru(ip1jmp1, llm), pbarv(ip1jm, llm)
105     REAL w(ip1jmp1, llm)
106     REAL phystep
107     INTEGER radpas
108    
109     !---------------------------------
110    
111     print *, "Call sequence information: etat0"
112    
113     ! Construct a grid:
114    
115     dtvr = daysec / real(day_step)
116     print *, 'dtvr = ', dtvr
117    
118     pa = 5e4
119     CALL iniconst
120     CALL inigeom
121     CALL inifilr
122    
123     latfi(1) = 90.
124     latfi(2:klon-1) = pack(spread(rlatu(2:jjm), 1, iim), .true.) * 180. / pi
125     ! (with conversion to degrees)
126     latfi(klon) = - 90.
127    
128     lonfi(1) = 0.
129     lonfi(2:klon-1) = pack(spread(rlonv(:iim), 2, jjm - 1), .true.) * 180. / pi
130     ! (with conversion to degrees)
131     lonfi(klon) = 0.
132    
133     call start_init_orog(relief, zstd_2d, zsig_2d, zgam_2d, zthe_2d, zpic_2d, &
134     zval_2d) ! also compute "masque" and "phis"
135     call init_dyn_phy ! define the mask "dyn_phy" for distinct grid points
136     zmasq = pack(masque, dyn_phy)
137     PRINT *, 'Masque construit'
138    
139     CALL start_init_dyn(tsol_2d, psol) ! also compute "qsol_2d"
140    
141     ! Compute pressure on intermediate levels:
142     forall(l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * psol(:, :)
143     CALL exner_hyb(psol, p3d, pks, pk)
144     IF (MINVAL(pk) == MAXVAL(pk)) stop '"pk" should not be constant'
145    
146     pls(:, :, :) = preff * (pk(:, :, :) / cpp)**(1. / kappa)
147     PRINT *, "minval(pls(:, :, :)) = ", minval(pls(:, :, :))
148     print *, "maxval(pls(:, :, :)) = ", maxval(pls(:, :, :))
149    
150     uvent(:, :, :) = start_inter_3d('U', rlonv, rlatv, pls)
151     forall (l = 1: llm) uvent(:iim, :, l) = uvent(:iim, :, l) * cu_2d(:iim, :)
152     uvent(iim+1, :, :) = uvent(1, :, :)
153    
154     vvent(:, :, :) = start_inter_3d('V', rlonu, rlatu(:jjm), pls(:, :jjm, :))
155     forall (l = 1: llm) vvent(:iim, :, l) = vvent(:iim, :, l) * cv_2d(:iim, :)
156     vvent(iim + 1, :, :) = vvent(1, :, :)
157    
158     t3d(:, :, :) = start_inter_3d('TEMP', rlonu, rlatv, pls)
159     PRINT *, 'minval(t3d(:, :, :)) = ', minval(t3d(:, :, :))
160     print *, "maxval(t3d(:, :, :)) = ", maxval(t3d(:, :, :))
161    
162     tpot(:iim, :, :) = t3d(:iim, :, :) * cpp / pk(:iim, :, :)
163     tpot(iim + 1, :, :) = tpot(1, :, :)
164     DO l=1, llm
165     tpot(:, 1, l) = SUM(aire_2d(:, 1) * tpot(:, 1, l)) / apoln
166     tpot(:, jjm + 1, l) = SUM(aire_2d(:, jjm + 1) * tpot(:, jjm + 1, l)) &
167     / apols
168     ENDDO
169    
170     ! Calcul de l'humidité à saturation :
171     qsat(:, :, :) = q_sat(t3d, pls)
172     PRINT *, "minval(qsat(:, :, :)) = ", minval(qsat(:, :, :))
173     print *, "maxval(qsat(:, :, :)) = ", maxval(qsat(:, :, :))
174     IF (MINVAL(qsat) == MAXVAL(qsat)) stop '"qsat" should not be constant'
175    
176     ! Water vapor:
177     q3d(:, :, :, 1) = 0.01 * start_inter_3d('R', rlonu, rlatv, pls) * qsat
178     WHERE (q3d(:, :, :, 1) < 0.) q3d(:, :, :, 1) = 1E-10
179     DO l = 1, llm
180     q3d(:, 1, l, 1) = SUM(aire_2d(:, 1) * q3d(:, 1, l, 1)) / apoln
181     q3d(:, jjm + 1, l, 1) &
182     = SUM(aire_2d(:, jjm + 1) * q3d(:, jjm + 1, l, 1)) / apols
183     ENDDO
184    
185     q3d(:, :, :, 2:4) = 0. ! liquid water, radon and lead
186    
187 guez 5 if (nqmx >= 5) then
188     ! Ozone:
189 guez 7 call regr_lat_time_coefoz
190     call regr_pr_o3(q3d(:, :, :, 5))
191     ! Convert from mole fraction to mass fraction:
192     q3d(:, :, :, 5) = q3d(:, :, :, 5) * 48. / 29.
193 guez 5 end if
194 guez 3
195     tsol(:) = pack(tsol_2d, dyn_phy)
196     qsol(:) = pack(qsol_2d, dyn_phy)
197     sn(:) = 0. ! snow
198     radsol(:) = 0.
199     tslab(:) = 0. ! IM "slab" ocean
200     seaice(:) = 0.
201     rugmer(:) = 0.001
202     zmea(:) = pack(relief, dyn_phy)
203     zstd(:) = pack(zstd_2d, dyn_phy)
204     zsig(:) = pack(zsig_2d, dyn_phy)
205     zgam(:) = pack(zgam_2d, dyn_phy)
206     zthe(:) = pack(zthe_2d, dyn_phy)
207     zpic(:) = pack(zpic_2d, dyn_phy)
208     zval(:) = pack(zval_2d, dyn_phy)
209    
210     rugsrel(:) = 0.
211     IF (ok_orodr) rugsrel(:) = MAX(1.e-05, zstd(:) * zsig(:) / 2)
212    
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     IF (MAXVAL( lon_lic(:, :) ) > pi) THEN
237     lon_lic(:, :) = lon_lic(:, :) * pi / 180.
238     ENDIF
239     IF (maxval( lat_lic(:, :) ) > pi) THEN
240     lat_lic(:, :) = lat_lic(:, :) * pi/ 180.
241     ENDIF
242    
243     dlon_lic(:) = lon_lic(:, 1)
244     dlat_lic(:) = lat_lic(1, :)
245    
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     pctsrf(:, :)=0.
252     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     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     pctsrf(:, is_ter) = 0.
261     elsewhere
262     pctsrf(:, is_ter) = zmasq(:) - pctsrf(:, is_lic)
263     where (pctsrf(:, is_ter) < EPSFRA)
264     pctsrf(:, is_ter) = 0.
265     pctsrf(:, is_lic) = zmasq(:)
266     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     pctsrf(:, is_oce) = 1. - zmasq(:)
273     WHERE (pctsrf(:, is_oce) < EPSFRA) pctsrf(:, is_oce) = 0.
274    
275     ! Vérification que somme des sous-surfaces vaut 1:
276     ji = count(abs(sum(pctsrf(:, :), dim = 2) - 1. ) > EPSFRA)
277     IF (ji /= 0) PRINT *, 'Problème répartition sous maille pour', ji, 'points'
278    
279     ! Calcul intermédiaire:
280     CALL massdair(p3d, masse)
281    
282     print *, 'ALPHAX = ', alphax
283    
284     forall (l = 1:llm)
285     masse(:, 1, l) = SUM(aire_2d(:iim, 1) * masse(:iim, 1, l)) / apoln
286     masse(:, jjm + 1, l) = &
287     SUM(aire_2d(:iim, jjm + 1) * masse(:iim, jjm + 1, l)) / apols
288     END forall
289    
290     ! Initialisation pour traceurs:
291 guez 5 call iniadvtrac
292 guez 3 ! Ecriture:
293     CALL inidissip(lstardis, nitergdiv, nitergrot, niterh, tetagdiv, &
294     tetagrot, tetatemp)
295     itau_dyn = 0
296     itau_phy = 0
297     day_ref = dayref
298     annee_ref = anneeref
299    
300     CALL geopot(ip1jmp1, tpot, pk , pks, phis , phi)
301     CALL caldyn0(0, uvent, vvent, tpot, psol, masse, pk, phis, phi, w, &
302     pbaru, pbarv, 0)
303 guez 5 CALL dynredem0("start.nc", dayref, phis)
304     CALL dynredem1("start.nc", 0., vvent, uvent, tpot, q3d, masse, psol)
305 guez 3
306     ! Ecriture état initial physique:
307     print *, 'dtvr = ', dtvr
308     print *, "iphysiq = ", iphysiq
309     print *, "nbapp_rad = ", nbapp_rad
310     phystep = dtvr * REAL(iphysiq)
311     radpas = NINT (86400./phystep/ nbapp_rad)
312     print *, 'phystep = ', phystep
313     print *, "radpas = ", radpas
314    
315     ! Initialisations :
316     tsolsrf(:, is_ter) = tsol
317     tsolsrf(:, is_lic) = tsol
318     tsolsrf(:, is_oce) = tsol
319     tsolsrf(:, is_sic) = tsol
320     snsrf(:, is_ter) = sn
321     snsrf(:, is_lic) = sn
322     snsrf(:, is_oce) = sn
323     snsrf(:, is_sic) = sn
324     albe(:, is_ter) = 0.08
325     albe(:, is_lic) = 0.6
326     albe(:, is_oce) = 0.5
327     albe(:, is_sic) = 0.6
328     alblw = albe
329     evap(:, :) = 0.
330     qsolsrf(:, is_ter) = 150.
331     qsolsrf(:, is_lic) = 150.
332     qsolsrf(:, is_oce) = 150.
333     qsolsrf(:, is_sic) = 150.
334     tsoil = spread(spread(tsol, 2, nsoilmx), 3, nbsrf)
335     rain_fall = 0.
336     snow_fall = 0.
337     solsw = 165.
338     sollw = -53.
339     t_ancien = 273.15
340     q_ancien = 0.
341     agesno = 0.
342     !IM "slab" ocean
343     tslab(:) = tsolsrf(:, is_oce)
344     seaice = 0.
345    
346     frugs(:, is_oce) = rugmer(:)
347     frugs(:, is_ter) = MAX(1.e-05, zstd(:) * zsig(:) / 2)
348     frugs(:, is_lic) = MAX(1.e-05, zstd(:) * zsig(:) / 2)
349     frugs(:, is_sic) = 0.001
350     fder = 0.
351     clwcon = 0.
352     rnebcon = 0.
353     ratqs = 0.
354     run_off_lic_0 = 0.
355    
356     call phyredem("startphy.nc", phystep, radpas, latfi, lonfi, pctsrf, &
357     tsolsrf, tsoil, tslab, seaice, qsolsrf, qsol, snsrf, albe, alblw, &
358     evap, rain_fall, snow_fall, solsw, sollw, fder, radsol, frugs, &
359     agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, rugsrel, &
360     t_ancien, q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0)
361     CALL histclo
362    
363     END SUBROUTINE etat0
364    
365     end module etat0_mod

  ViewVC Help
Powered by ViewVC 1.1.21