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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5 - (hide annotations)
Mon Mar 3 16:32:04 2008 UTC (16 years, 2 months ago) by guez
File size: 13035 byte(s)
Created module from included file parafilt.
Converted caldyn0 to free format.
Added a rule to create cross-references with NAG.
Added optional attribute in iniadvtrac.
Suppressed argument nq in dynredem0 and dynredem1, using nqmx instead.
Replaced some NetCDF calls by netcdf95 calls in dynredem0.
Added intent attribute in dynredem0 and dynredem1.
Annotated use statements with only clause, in dynredem1.
Suppressed variable nq and argument of iniadvtrac in etat0.
Added test on nqmx in etat0.
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    
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     integer ncid, varid, ncerr, month
109    
110     !---------------------------------
111    
112     print *, "Call sequence information: etat0"
113    
114     ! Construct a grid:
115    
116     dtvr = daysec / real(day_step)
117     print *, 'dtvr = ', dtvr
118    
119     pa = 5e4
120     CALL iniconst
121     CALL inigeom
122     CALL inifilr
123    
124     latfi(1) = 90.
125     latfi(2:klon-1) = pack(spread(rlatu(2:jjm), 1, iim), .true.) * 180. / pi
126     ! (with conversion to degrees)
127     latfi(klon) = - 90.
128    
129     lonfi(1) = 0.
130     lonfi(2:klon-1) = pack(spread(rlonv(:iim), 2, jjm - 1), .true.) * 180. / pi
131     ! (with conversion to degrees)
132     lonfi(klon) = 0.
133    
134     call start_init_orog(relief, zstd_2d, zsig_2d, zgam_2d, zthe_2d, zpic_2d, &
135     zval_2d) ! also compute "masque" and "phis"
136     call init_dyn_phy ! define the mask "dyn_phy" for distinct grid points
137     zmasq = pack(masque, dyn_phy)
138     PRINT *, 'Masque construit'
139    
140     CALL start_init_dyn(tsol_2d, psol) ! also compute "qsol_2d"
141    
142     ! Compute pressure on intermediate levels:
143     forall(l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * psol(:, :)
144     CALL exner_hyb(psol, p3d, pks, pk)
145     IF (MINVAL(pk) == MAXVAL(pk)) stop '"pk" should not be constant'
146    
147     pls(:, :, :) = preff * (pk(:, :, :) / cpp)**(1. / kappa)
148     PRINT *, "minval(pls(:, :, :)) = ", minval(pls(:, :, :))
149     print *, "maxval(pls(:, :, :)) = ", maxval(pls(:, :, :))
150    
151     uvent(:, :, :) = start_inter_3d('U', rlonv, rlatv, pls)
152     forall (l = 1: llm) uvent(:iim, :, l) = uvent(:iim, :, l) * cu_2d(:iim, :)
153     uvent(iim+1, :, :) = uvent(1, :, :)
154    
155     vvent(:, :, :) = start_inter_3d('V', rlonu, rlatu(:jjm), pls(:, :jjm, :))
156     forall (l = 1: llm) vvent(:iim, :, l) = vvent(:iim, :, l) * cv_2d(:iim, :)
157     vvent(iim + 1, :, :) = vvent(1, :, :)
158    
159     t3d(:, :, :) = start_inter_3d('TEMP', rlonu, rlatv, pls)
160     PRINT *, 'minval(t3d(:, :, :)) = ', minval(t3d(:, :, :))
161     print *, "maxval(t3d(:, :, :)) = ", maxval(t3d(:, :, :))
162    
163     tpot(:iim, :, :) = t3d(:iim, :, :) * cpp / pk(:iim, :, :)
164     tpot(iim + 1, :, :) = tpot(1, :, :)
165     DO l=1, llm
166     tpot(:, 1, l) = SUM(aire_2d(:, 1) * tpot(:, 1, l)) / apoln
167     tpot(:, jjm + 1, l) = SUM(aire_2d(:, jjm + 1) * tpot(:, jjm + 1, l)) &
168     / apols
169     ENDDO
170    
171     ! Calcul de l'humidité à saturation :
172     qsat(:, :, :) = q_sat(t3d, pls)
173     PRINT *, "minval(qsat(:, :, :)) = ", minval(qsat(:, :, :))
174     print *, "maxval(qsat(:, :, :)) = ", maxval(qsat(:, :, :))
175     IF (MINVAL(qsat) == MAXVAL(qsat)) stop '"qsat" should not be constant'
176    
177     ! Water vapor:
178     q3d(:, :, :, 1) = 0.01 * start_inter_3d('R', rlonu, rlatv, pls) * qsat
179     WHERE (q3d(:, :, :, 1) < 0.) q3d(:, :, :, 1) = 1E-10
180     DO l = 1, llm
181     q3d(:, 1, l, 1) = SUM(aire_2d(:, 1) * q3d(:, 1, l, 1)) / apoln
182     q3d(:, jjm + 1, l, 1) &
183     = SUM(aire_2d(:, jjm + 1) * q3d(:, jjm + 1, l, 1)) / apols
184     ENDDO
185    
186     q3d(:, :, :, 2:4) = 0. ! liquid water, radon and lead
187    
188 guez 5 if (nqmx >= 5) then
189     ! Ozone:
190 guez 3
191 guez 5 ! Compute ozone parameters on the LMDZ grid:
192     call regr_coefoz
193 guez 3
194 guez 5 ! Find the month containing the day number "dayref":
195     month = (dayref - 1) / 30 + 1
196     print *, "month = ", month
197 guez 3
198 guez 5 call nf95_open("coefoz_LMDZ.nc", nf90_nowrite, ncid)
199 guez 3
200 guez 5 ! 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 guez 3
206 guez 5 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     end if
212 guez 3
213     tsol(:) = pack(tsol_2d, dyn_phy)
214     qsol(:) = pack(qsol_2d, dyn_phy)
215     sn(:) = 0. ! snow
216     radsol(:) = 0.
217     tslab(:) = 0. ! IM "slab" ocean
218     seaice(:) = 0.
219     rugmer(:) = 0.001
220     zmea(:) = pack(relief, dyn_phy)
221     zstd(:) = pack(zstd_2d, dyn_phy)
222     zsig(:) = pack(zsig_2d, dyn_phy)
223     zgam(:) = pack(zgam_2d, dyn_phy)
224     zthe(:) = pack(zthe_2d, dyn_phy)
225     zpic(:) = pack(zpic_2d, dyn_phy)
226     zval(:) = pack(zval_2d, dyn_phy)
227    
228     rugsrel(:) = 0.
229     IF (ok_orodr) rugsrel(:) = MAX(1.e-05, zstd(:) * zsig(:) / 2)
230    
231     ! On initialise les sous-surfaces:
232     ! Lecture du fichier glace de terre pour fixer la fraction de terre
233     ! et de glace de terre:
234     CALL flininfo("landiceref.nc", iml_lic, jml_lic, llm_tmp, &
235     ttm_tmp, fid)
236     ALLOCATE(lat_lic(iml_lic, jml_lic))
237     ALLOCATE(lon_lic(iml_lic, jml_lic))
238     ALLOCATE(dlon_lic(iml_lic))
239     ALLOCATE(dlat_lic(jml_lic))
240     ALLOCATE(fraclic(iml_lic, jml_lic))
241     CALL flinopen_nozoom("landiceref.nc", iml_lic, jml_lic, &
242     llm_tmp, lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, &
243     fid)
244     CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp &
245     , 1, 1, fraclic)
246     CALL flinclo(fid)
247    
248     ! Interpolation sur la grille T du modèle :
249     PRINT *, 'Dimensions de "landice"'
250     print *, "iml_lic = ", iml_lic
251     print *, "jml_lic = ", jml_lic
252    
253     ! Si les coordonnées sont en degrés, on les transforme :
254     IF (MAXVAL( lon_lic(:, :) ) > pi) THEN
255     lon_lic(:, :) = lon_lic(:, :) * pi / 180.
256     ENDIF
257     IF (maxval( lat_lic(:, :) ) > pi) THEN
258     lat_lic(:, :) = lat_lic(:, :) * pi/ 180.
259     ENDIF
260    
261     dlon_lic(:) = lon_lic(:, 1)
262     dlat_lic(:) = lat_lic(1, :)
263    
264     flic_tmp(:iim, :) = grille_m(dlon_lic, dlat_lic, fraclic, rlonv(:iim), &
265     rlatu)
266     flic_tmp(iim + 1, :) = flic_tmp(1, :)
267    
268     ! Passage sur la grille physique
269     pctsrf(:, :)=0.
270     pctsrf(:, is_lic) = pack(flic_tmp, dyn_phy)
271     ! Adéquation avec le maque terre/mer
272     WHERE (pctsrf(:, is_lic) < EPSFRA ) pctsrf(:, is_lic) = 0.
273     WHERE (zmasq(:) < EPSFRA) pctsrf(:, is_lic) = 0.
274     pctsrf(:, is_ter) = zmasq(:)
275     where (zmasq(:) > EPSFRA)
276     where (pctsrf(:, is_lic) >= zmasq(:))
277     pctsrf(:, is_lic) = zmasq(:)
278     pctsrf(:, is_ter) = 0.
279     elsewhere
280     pctsrf(:, is_ter) = zmasq(:) - pctsrf(:, is_lic)
281     where (pctsrf(:, is_ter) < EPSFRA)
282     pctsrf(:, is_ter) = 0.
283     pctsrf(:, is_lic) = zmasq(:)
284     end where
285     end where
286     end where
287    
288     ! Sous-surface océan et glace de mer (pour démarrer on met glace
289     ! de mer à 0) :
290     pctsrf(:, is_oce) = 1. - zmasq(:)
291     WHERE (pctsrf(:, is_oce) < EPSFRA) pctsrf(:, is_oce) = 0.
292    
293     ! Vérification que somme des sous-surfaces vaut 1:
294     ji = count(abs(sum(pctsrf(:, :), dim = 2) - 1. ) > EPSFRA)
295     IF (ji /= 0) PRINT *, 'Problème répartition sous maille pour', ji, 'points'
296    
297     ! Calcul intermédiaire:
298     CALL massdair(p3d, masse)
299    
300     print *, 'ALPHAX = ', alphax
301    
302     forall (l = 1:llm)
303     masse(:, 1, l) = SUM(aire_2d(:iim, 1) * masse(:iim, 1, l)) / apoln
304     masse(:, jjm + 1, l) = &
305     SUM(aire_2d(:iim, jjm + 1) * masse(:iim, jjm + 1, l)) / apols
306     END forall
307    
308     ! Initialisation pour traceurs:
309 guez 5 call iniadvtrac
310 guez 3 ! Ecriture:
311     CALL inidissip(lstardis, nitergdiv, nitergrot, niterh, tetagdiv, &
312     tetagrot, tetatemp)
313     itau_dyn = 0
314     itau_phy = 0
315     day_ref = dayref
316     annee_ref = anneeref
317    
318     CALL geopot(ip1jmp1, tpot, pk , pks, phis , phi)
319     CALL caldyn0(0, uvent, vvent, tpot, psol, masse, pk, phis, phi, w, &
320     pbaru, pbarv, 0)
321 guez 5 CALL dynredem0("start.nc", dayref, phis)
322     CALL dynredem1("start.nc", 0., vvent, uvent, tpot, q3d, masse, psol)
323 guez 3
324     ! Ecriture état initial physique:
325     print *, 'dtvr = ', dtvr
326     print *, "iphysiq = ", iphysiq
327     print *, "nbapp_rad = ", nbapp_rad
328     phystep = dtvr * REAL(iphysiq)
329     radpas = NINT (86400./phystep/ nbapp_rad)
330     print *, 'phystep = ', phystep
331     print *, "radpas = ", radpas
332    
333     ! Initialisations :
334     tsolsrf(:, is_ter) = tsol
335     tsolsrf(:, is_lic) = tsol
336     tsolsrf(:, is_oce) = tsol
337     tsolsrf(:, is_sic) = tsol
338     snsrf(:, is_ter) = sn
339     snsrf(:, is_lic) = sn
340     snsrf(:, is_oce) = sn
341     snsrf(:, is_sic) = sn
342     albe(:, is_ter) = 0.08
343     albe(:, is_lic) = 0.6
344     albe(:, is_oce) = 0.5
345     albe(:, is_sic) = 0.6
346     alblw = albe
347     evap(:, :) = 0.
348     qsolsrf(:, is_ter) = 150.
349     qsolsrf(:, is_lic) = 150.
350     qsolsrf(:, is_oce) = 150.
351     qsolsrf(:, is_sic) = 150.
352     tsoil = spread(spread(tsol, 2, nsoilmx), 3, nbsrf)
353     rain_fall = 0.
354     snow_fall = 0.
355     solsw = 165.
356     sollw = -53.
357     t_ancien = 273.15
358     q_ancien = 0.
359     agesno = 0.
360     !IM "slab" ocean
361     tslab(:) = tsolsrf(:, is_oce)
362     seaice = 0.
363    
364     frugs(:, is_oce) = rugmer(:)
365     frugs(:, is_ter) = MAX(1.e-05, zstd(:) * zsig(:) / 2)
366     frugs(:, is_lic) = MAX(1.e-05, zstd(:) * zsig(:) / 2)
367     frugs(:, is_sic) = 0.001
368     fder = 0.
369     clwcon = 0.
370     rnebcon = 0.
371     ratqs = 0.
372     run_off_lic_0 = 0.
373    
374     call phyredem("startphy.nc", phystep, radpas, latfi, lonfi, pctsrf, &
375     tsolsrf, tsoil, tslab, seaice, qsolsrf, qsol, snsrf, albe, alblw, &
376     evap, rain_fall, snow_fall, solsw, sollw, fder, radsol, frugs, &
377     agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, rugsrel, &
378     t_ancien, q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0)
379     CALL histclo
380    
381     END SUBROUTINE etat0
382    
383     end module etat0_mod

  ViewVC Help
Powered by ViewVC 1.1.21