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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 18 - (hide annotations)
Thu Aug 7 12:29:13 2008 UTC (15 years, 9 months ago) by guez
File size: 12017 byte(s)
In module "regr_pr", rewrote scanning of horizontal positions as a
single set of loops, using a mask.

Added some "intent" attributes.

In "dynredem0", replaced calls to Fortran 77 interface of NetCDF by
calls to NetCDF95. Removed calls to "nf_redef", regrouped all writing
operations. In "dynredem1", replaced some calls to Fortran 77
interface of NetCDF by calls to Fortran 90 interface.

Renamed variable "nqmax" to "nq_phys".

In "physiq", if "nq >= 5" then "wo" is computed from the
parameterization of "Cariolle".

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 3
52     ! Variables local to the procedure:
53    
54     REAL latfi(klon), lonfi(klon)
55     ! (latitude and longitude of a point of the scalar grid identified
56     ! by a simple index, in °)
57    
58     REAL, dimension(iim + 1, jjm + 1, llm):: uvent, t3d, tpot
59     REAL vvent(iim + 1, jjm, llm)
60    
61     REAL q3d(iim + 1, jjm + 1, llm, nqmx)
62     ! (mass fractions of trace species
63     ! "q3d(i, j, l)" is at longitude "rlonv(i)", latitude "rlatu(j)"
64     ! and pressure level "pls(i, j, l)".)
65    
66     real qsat(iim + 1, jjm + 1, llm) ! mass fraction of saturating water vapor
67     REAL tsol(klon), qsol(klon), sn(klon)
68     REAL tsolsrf(klon, nbsrf), qsolsrf(klon, nbsrf), snsrf(klon, nbsrf)
69     REAL albe(klon, nbsrf), evap(klon, nbsrf)
70     REAL alblw(klon, nbsrf)
71     REAL tsoil(klon, nsoilmx, nbsrf)
72     REAL radsol(klon), rain_fall(klon), snow_fall(klon)
73     REAL solsw(klon), sollw(klon), fder(klon)
74     !IM "slab" ocean
75     REAL tslab(klon)
76     real seaice(klon) ! kg m-2
77     REAL frugs(klon, nbsrf), agesno(klon, nbsrf)
78     REAL rugmer(klon)
79     real, dimension(iim + 1, jjm + 1):: relief, zstd_2d, zsig_2d, zgam_2d
80     real, dimension(iim + 1, jjm + 1):: zthe_2d, zpic_2d, zval_2d
81     real, dimension(iim + 1, jjm + 1):: tsol_2d, psol
82     REAL zmea(klon), zstd(klon)
83     REAL zsig(klon), zgam(klon)
84     REAL zthe(klon)
85     REAL zpic(klon), zval(klon)
86     REAL t_ancien(klon, llm), q_ancien(klon, llm) !
87     REAL run_off_lic_0(klon)
88     real clwcon(klon, llm), rnebcon(klon, llm), ratqs(klon, llm)
89     ! déclarations pour lecture glace de mer
90     INTEGER iml_lic, jml_lic, llm_tmp, ttm_tmp
91     INTEGER itaul(1), fid
92     REAL lev(1), date
93     REAL, ALLOCATABLE:: lon_lic(:, :), lat_lic(:, :)
94     REAL, ALLOCATABLE:: dlon_lic(:), dlat_lic(:)
95     REAL, ALLOCATABLE:: fraclic(:, :) ! fraction land ice
96     REAL flic_tmp(iim + 1, jjm + 1) !fraction land ice temporary
97    
98     INTEGER l, ji
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    
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 guez 15 zval_2d) ! also compute "mask" and "phis"
135 guez 3 call init_dyn_phy ! define the mask "dyn_phy" for distinct grid points
136 guez 15 zmasq = pack(mask, dyn_phy)
137 guez 3 PRINT *, 'Masque construit'
138    
139     CALL start_init_dyn(tsol_2d, psol) ! also compute "qsol_2d"
140    
141     ! Compute pressure on intermediate levels:
142 guez 15 forall(l = 1: llm + 1) p3d(:, :, l) = ap(l) + bp(l) * psol
143 guez 3 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 guez 13 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 guez 3
210     ! On initialise les sous-surfaces:
211     ! Lecture du fichier glace de terre pour fixer la fraction de terre
212     ! et de glace de terre:
213     CALL flininfo("landiceref.nc", iml_lic, jml_lic, llm_tmp, &
214     ttm_tmp, fid)
215     ALLOCATE(lat_lic(iml_lic, jml_lic))
216     ALLOCATE(lon_lic(iml_lic, jml_lic))
217     ALLOCATE(dlon_lic(iml_lic))
218     ALLOCATE(dlat_lic(jml_lic))
219     ALLOCATE(fraclic(iml_lic, jml_lic))
220     CALL flinopen_nozoom("landiceref.nc", iml_lic, jml_lic, &
221     llm_tmp, lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, &
222     fid)
223     CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp &
224     , 1, 1, fraclic)
225     CALL flinclo(fid)
226    
227     ! Interpolation sur la grille T du modèle :
228     PRINT *, 'Dimensions de "landice"'
229     print *, "iml_lic = ", iml_lic
230     print *, "jml_lic = ", jml_lic
231    
232     ! Si les coordonnées sont en degrés, on les transforme :
233 guez 15 IF (MAXVAL( lon_lic ) > pi) THEN
234     lon_lic = lon_lic * pi / 180.
235 guez 3 ENDIF
236 guez 15 IF (maxval( lat_lic ) > pi) THEN
237     lat_lic = lat_lic * pi/ 180.
238 guez 3 ENDIF
239    
240 guez 13 dlon_lic = lon_lic(:, 1)
241     dlat_lic = lat_lic(1, :)
242 guez 3
243     flic_tmp(:iim, :) = grille_m(dlon_lic, dlat_lic, fraclic, rlonv(:iim), &
244     rlatu)
245     flic_tmp(iim + 1, :) = flic_tmp(1, :)
246    
247     ! Passage sur la grille physique
248 guez 15 pctsrf = 0.
249 guez 3 pctsrf(:, is_lic) = pack(flic_tmp, dyn_phy)
250     ! Adéquation avec le maque terre/mer
251     WHERE (pctsrf(:, is_lic) < EPSFRA ) pctsrf(:, is_lic) = 0.
252 guez 13 WHERE (zmasq < EPSFRA) pctsrf(:, is_lic) = 0.
253     pctsrf(:, is_ter) = zmasq
254     where (zmasq > EPSFRA)
255     where (pctsrf(:, is_lic) >= zmasq)
256     pctsrf(:, is_lic) = zmasq
257 guez 3 pctsrf(:, is_ter) = 0.
258     elsewhere
259 guez 13 pctsrf(:, is_ter) = zmasq - pctsrf(:, is_lic)
260 guez 3 where (pctsrf(:, is_ter) < EPSFRA)
261     pctsrf(:, is_ter) = 0.
262 guez 13 pctsrf(:, is_lic) = zmasq
263 guez 3 end where
264     end where
265     end where
266    
267     ! Sous-surface océan et glace de mer (pour démarrer on met glace
268     ! de mer à 0) :
269 guez 13 pctsrf(:, is_oce) = 1. - zmasq
270 guez 3 WHERE (pctsrf(:, is_oce) < EPSFRA) pctsrf(:, is_oce) = 0.
271    
272     ! Vérification que somme des sous-surfaces vaut 1:
273 guez 15 ji = count(abs(sum(pctsrf, dim = 2) - 1.) > EPSFRA)
274     IF (ji /= 0) then
275     PRINT *, 'Problème répartition sous maille pour ', ji, 'points'
276     end IF
277 guez 3
278     ! Calcul intermédiaire:
279     CALL massdair(p3d, masse)
280    
281     print *, 'ALPHAX = ', alphax
282    
283     forall (l = 1:llm)
284     masse(:, 1, l) = SUM(aire_2d(:iim, 1) * masse(:iim, 1, l)) / apoln
285     masse(:, jjm + 1, l) = &
286     SUM(aire_2d(:iim, jjm + 1) * masse(:iim, jjm + 1, l)) / apols
287     END forall
288    
289     ! Initialisation pour traceurs:
290 guez 5 call iniadvtrac
291 guez 3 ! Ecriture:
292     CALL inidissip(lstardis, nitergdiv, nitergrot, niterh, tetagdiv, &
293     tetagrot, tetatemp)
294     itau_dyn = 0
295     itau_phy = 0
296     day_ref = dayref
297     annee_ref = anneeref
298    
299     CALL geopot(ip1jmp1, tpot, pk , pks, phis , phi)
300     CALL caldyn0(0, uvent, vvent, tpot, psol, masse, pk, phis, phi, w, &
301     pbaru, pbarv, 0)
302 guez 5 CALL dynredem0("start.nc", dayref, phis)
303     CALL dynredem1("start.nc", 0., vvent, uvent, tpot, q3d, masse, psol)
304 guez 3
305     ! Ecriture état initial physique:
306     print *, 'dtvr = ', dtvr
307     print *, "iphysiq = ", iphysiq
308     phystep = dtvr * REAL(iphysiq)
309     print *, 'phystep = ', phystep
310    
311     ! Initialisations :
312     tsolsrf(:, is_ter) = tsol
313     tsolsrf(:, is_lic) = tsol
314     tsolsrf(:, is_oce) = tsol
315     tsolsrf(:, is_sic) = tsol
316     snsrf(:, is_ter) = sn
317     snsrf(:, is_lic) = sn
318     snsrf(:, is_oce) = sn
319     snsrf(:, is_sic) = sn
320     albe(:, is_ter) = 0.08
321     albe(:, is_lic) = 0.6
322     albe(:, is_oce) = 0.5
323     albe(:, is_sic) = 0.6
324     alblw = albe
325 guez 15 evap = 0.
326 guez 3 qsolsrf(:, is_ter) = 150.
327     qsolsrf(:, is_lic) = 150.
328     qsolsrf(:, is_oce) = 150.
329     qsolsrf(:, is_sic) = 150.
330     tsoil = spread(spread(tsol, 2, nsoilmx), 3, nbsrf)
331     rain_fall = 0.
332     snow_fall = 0.
333     solsw = 165.
334     sollw = -53.
335     t_ancien = 273.15
336     q_ancien = 0.
337     agesno = 0.
338     !IM "slab" ocean
339 guez 13 tslab = tsolsrf(:, is_oce)
340 guez 3 seaice = 0.
341    
342 guez 13 frugs(:, is_oce) = rugmer
343     frugs(:, is_ter) = MAX(1.e-05, zstd * zsig / 2)
344     frugs(:, is_lic) = MAX(1.e-05, zstd * zsig / 2)
345 guez 3 frugs(:, is_sic) = 0.001
346     fder = 0.
347     clwcon = 0.
348     rnebcon = 0.
349     ratqs = 0.
350     run_off_lic_0 = 0.
351    
352 guez 13 call phyredem("startphy.nc", latfi, lonfi, pctsrf, &
353 guez 3 tsolsrf, tsoil, tslab, seaice, qsolsrf, qsol, snsrf, albe, alblw, &
354     evap, rain_fall, snow_fall, solsw, sollw, fder, radsol, frugs, &
355 guez 13 agesno, zmea, zstd, zsig, zgam, zthe, zpic, zval, &
356 guez 3 t_ancien, q_ancien, rnebcon, ratqs, clwcon, run_off_lic_0)
357     CALL histclo
358    
359     END SUBROUTINE etat0
360    
361     end module etat0_mod

  ViewVC Help
Powered by ViewVC 1.1.21