/[lmdze]/trunk/dyn3d/limit.f
ViewVC logotype

Contents of /trunk/dyn3d/limit.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 98 - (show annotations)
Tue May 13 17:23:16 2014 UTC (10 years ago) by guez
File size: 13909 byte(s)
Split inter_barxy.f : one procedure per module, one module per
file. Grouped the files into a directory.

Split orbite.f.

Value of raz_date read from the namelist is taken into account
(resetting the step counter) even if annee_ref == anneeref and day_ref
== dayref. raz_date is no longer modified by gcm main unit. (Following
LMDZ.)

Removed argument klon of interfsur_lim. Renamed arguments lmt_alb,
lmt_rug to alb_new, z0_new (same name as corresponding actual
arguments in interfsurf_hq).

Removed argument klon of interfsurf_hq.

Removed arguments qs and d_qs of diagetpq. Were always
zero. Downgraded arguments d_qw, d_ql of diagetpq to local variables,
they were not used in physiq. Removed all computations for solid water
in diagetpq, was just zero.


Downgraded arguments fs_bound, fq_bound of diagphy to local variables,
they were not used in physiq. Encapsulated in a test on iprt all
computations in diagphy.

Removed parameter nbtr of module dimphy. Replaced it everywhere in the
program by nqmx - 2.

Removed parameter rnpb of procedure physiq. Kept the true case in
physiq and phytrac. Could not work with false case anyway.

Removed arguments klon, llm, airephy of qcheck. Removed argument ftsol
of initrrnpb, was not used.

1 module limit_mod
2
3 IMPLICIT none
4
5 contains
6
7 SUBROUTINE limit
8
9 ! Authors: L. Fairhead, Z. X. Li, P. Le Van
10
11 ! This subroutine creates files containing boundary conditions.
12 ! It uses files with climatological data.
13 ! Both grids must be regular.
14
15 use comgeom, only: rlonu, rlatv
16 use conf_dat2d_m, only: conf_dat2d
17 use dimens_m, only: iim, jjm
18 use dimphy, only: klon, zmasq
19 use etat0_mod, only: pctsrf
20 use grid_change, only: dyn_phy
21 use indicesol, only: epsfra, nbsrf, is_ter, is_oce, is_lic, is_sic
22 use inter_barxy_m, only: inter_barxy
23 use netcdf95, only: NF95_CLOSE, NF95_CREATE, NF95_DEF_DIM, nf95_def_var, &
24 nf95_enddef, nf95_get_var, nf95_gw_var, nf95_inq_dimid, &
25 nf95_inq_varid, nf95_inquire_dimension, NF95_OPEN, NF95_PUT_ATT, &
26 NF95_PUT_VAR
27 use netcdf, only: NF90_CLOBBER, NF90_FLOAT, NF90_GLOBAL, NF90_NOWRITE, &
28 NF90_UNLIMITED
29 use numer_rec_95, only: spline, splint
30 use start_init_orog_m, only: mask
31 use unit_nml_m, only: unit_nml
32
33 ! Variables local to the procedure:
34
35 LOGICAL:: extrap = .FALSE.
36 ! (extrapolation de données, comme pour les SST lorsque le fichier
37 ! ne contient pas uniquement des points océaniques)
38
39 REAL phy_alb(klon, 360)
40 REAL phy_sst(klon, 360)
41 REAL phy_bil(klon, 360)
42 REAL phy_rug(klon, 360)
43 REAL phy_ice(klon)
44
45 real pctsrf_t(klon, nbsrf, 360) ! composition of the surface
46
47 ! Pour le champ de départ:
48 INTEGER imdep, jmdep, lmdep
49
50 REAL, ALLOCATABLE:: dlon(:), dlat(:)
51 REAL, pointer:: dlon_ini(:), dlat_ini(:), timeyear(:)
52 REAL, ALLOCATABLE:: champ(:, :)
53 REAL, ALLOCATABLE:: work(:, :)
54
55 ! Pour le champ interpolé 3D :
56 REAL, allocatable:: champtime(:, :, :)
57 REAL champan(iim + 1, jjm + 1, 360)
58
59 ! Pour l'inteprolation verticale :
60 REAL, allocatable:: yder(:)
61
62 INTEGER nid, ndim, ntim
63 INTEGER id_tim
64 INTEGER id_SST, id_BILS, id_RUG, id_ALB
65 INTEGER id_FOCE, id_FSIC, id_FTER, id_FLIC
66
67 INTEGER i, j, k, l
68 INTEGER ncid, varid, dimid
69
70 REAL, parameter:: tmidmonth(12) = (/(15. + 30. * i, i = 0, 11)/)
71
72 namelist /limit_nml/extrap
73
74 !--------------------
75
76 print *, "Call sequence information: limit"
77
78 print *, "Enter namelist 'limit_nml'."
79 read(unit=*, nml=limit_nml)
80 write(unit_nml, nml=limit_nml)
81
82 PRINT *, 'Processing rugosity...'
83
84 call NF95_OPEN('Rugos.nc', NF90_NOWRITE, ncid)
85
86 ! Read coordinate variables:
87
88 call nf95_inq_varid(ncid, "longitude", varid)
89 call nf95_gw_var(ncid, varid, dlon_ini)
90 imdep = size(dlon_ini)
91
92 call nf95_inq_varid(ncid, "latitude", varid)
93 call nf95_gw_var(ncid, varid, dlat_ini)
94 jmdep = size(dlat_ini)
95
96 call nf95_inq_varid(ncid, "temps", varid)
97 call nf95_gw_var(ncid, varid, timeyear)
98 lmdep = size(timeyear)
99
100 ALLOCATE(champ(imdep, jmdep), champtime(iim, jjm + 1, lmdep))
101 allocate(dlon(imdep), dlat(jmdep))
102 call NF95_INQ_VARID(ncid, 'RUGOS', varid)
103
104 ! Read the primary variable day by day and regrid horizontally,
105 ! result in "champtime":
106 DO l = 1, lmdep
107 call NF95_GET_VAR(ncid, varid, champ, start=(/1, 1, l/))
108 CALL conf_dat2d(dlon_ini, dlat_ini, dlon, dlat, champ)
109 CALL inter_barxy(dlon, dlat(:jmdep -1), LOG(champ), rlonu(:iim), &
110 rlatv, champtime(:, :, l))
111 champtime(:, :, l) = EXP(champtime(:, :, l))
112 where (nint(mask(:iim, :)) /= 1) champtime(:, :, l) = 0.001
113 end do
114
115 call NF95_CLOSE(ncid)
116
117 DEALLOCATE(dlon, dlat, champ, dlon_ini, dlat_ini)
118 allocate(yder(lmdep))
119
120 ! Interpolate monthly values to daily values, at each horizontal position:
121 DO j = 1, jjm + 1
122 DO i = 1, iim
123 yder = SPLINE(timeyear, champtime(i, j, :))
124 DO k = 1, 360
125 champan(i, j, k) = SPLINT(timeyear, champtime(i, j, :), yder, &
126 real(k-1))
127 ENDDO
128 ENDDO
129 ENDDO
130
131 deallocate(timeyear, champtime, yder)
132 champan(iim + 1, :, :) = champan(1, :, :)
133 forall (k = 1:360) phy_rug(:, k) = pack(champan(:, :, k), dyn_phy)
134
135 ! Process sea ice:
136
137 PRINT *, 'Processing sea ice...'
138 call NF95_OPEN('amipbc_sic_1x1.nc', NF90_NOWRITE, ncid)
139
140 call nf95_inq_varid(ncid, "longitude", varid)
141 call nf95_gw_var(ncid, varid, dlon_ini)
142 imdep = size(dlon_ini)
143
144 call nf95_inq_varid(ncid, "latitude", varid)
145 call nf95_gw_var(ncid, varid, dlat_ini)
146 jmdep = size(dlat_ini)
147
148 call nf95_inq_dimid(ncid, "time", dimid)
149 call NF95_INQuire_DIMension(ncid, dimid, nclen=lmdep)
150 print *, 'lmdep = ', lmdep
151 ! Coordonnée temporelle fichiers AMIP pas en jours. Ici on suppose
152 ! qu'on a 12 mois (de 30 jours).
153 IF (lmdep /= 12) then
154 print *, 'Unknown AMIP file: not 12 months?'
155 STOP 1
156 end IF
157
158 ALLOCATE(champ(imdep, jmdep), champtime(iim, jjm + 1, lmdep))
159 ALLOCATE(dlon(imdep), dlat(jmdep))
160 call NF95_INQ_VARID(ncid, 'sicbcs', varid)
161 DO l = 1, lmdep
162 call NF95_GET_VAR(ncid, varid, champ, start=(/1, 1, l/))
163 CALL conf_dat2d(dlon_ini, dlat_ini, dlon, dlat, champ)
164 CALL inter_barxy(dlon, dlat(:jmdep -1), champ, rlonu(:iim), rlatv, &
165 champtime(:, :, l))
166 ENDDO
167
168 call NF95_CLOSE(ncid)
169
170 DEALLOCATE(dlon, dlat, champ, dlon_ini, dlat_ini)
171 PRINT *, 'Interpolation temporelle'
172 allocate(yder(lmdep))
173
174 DO j = 1, jjm + 1
175 DO i = 1, iim
176 yder = SPLINE(tmidmonth, champtime(i, j, :))
177 DO k = 1, 360
178 champan(i, j, k) = SPLINT(tmidmonth, champtime(i, j, :), yder, &
179 real(k-1))
180 ENDDO
181 ENDDO
182 ENDDO
183
184 deallocate(champtime, yder)
185
186 ! Convert from percentage to normal fraction and keep sea ice
187 ! between 0 and 1:
188 champan(:iim, :, :) = max(0., (min(1., (champan(:iim, :, :) / 100.))))
189 champan(iim + 1, :, :) = champan(1, :, :)
190
191 DO k = 1, 360
192 phy_ice = pack(champan(:, :, k), dyn_phy)
193
194 ! (utilisation de la sous-maille fractionnelle tandis que l'ancien
195 ! codage utilisait l'indicateur du sol (0, 1, 2, 3))
196 ! PB en attendant de mettre fraction de terre
197 WHERE (phy_ice < EPSFRA) phy_ice = 0.
198
199 pctsrf_t(:, is_ter, k) = pctsrf(:, is_ter)
200 pctsrf_t(:, is_lic, k) = pctsrf(:, is_lic)
201 pctsrf_t(:, is_sic, k) = max(phy_ice - pctsrf_t(:, is_lic, k), 0.)
202 ! Il y a des cas où il y a de la glace dans landiceref et
203 ! pas dans AMIP
204 WHERE (1. - zmasq < EPSFRA)
205 pctsrf_t(:, is_sic, k) = 0.
206 pctsrf_t(:, is_oce, k) = 0.
207 elsewhere
208 where (pctsrf_t(:, is_sic, k) >= 1 - zmasq)
209 pctsrf_t(:, is_sic, k) = 1. - zmasq
210 pctsrf_t(:, is_oce, k) = 0.
211 ELSEwhere
212 pctsrf_t(:, is_oce, k) = 1. - zmasq - pctsrf_t(:, is_sic, k)
213 where (pctsrf_t(:, is_oce, k) < EPSFRA)
214 pctsrf_t(:, is_oce, k) = 0.
215 pctsrf_t(:, is_sic, k) = 1 - zmasq
216 end where
217 end where
218 end where
219
220 DO i = 1, klon
221 if (pctsrf_t(i, is_oce, k) < 0.) then
222 print *, 'Problème sous maille : pctsrf_t(', i, &
223 ', is_oce, ', k, ') = ', pctsrf_t(i, is_oce, k)
224 ENDIF
225 IF (abs(pctsrf_t(i, is_ter, k) + pctsrf_t(i, is_lic, k) &
226 + pctsrf_t(i, is_oce, k) + pctsrf_t(i, is_sic, k) - 1.) &
227 > EPSFRA) THEN
228 print *, 'Problème sous surface :'
229 print *, "pctsrf_t(", i, ", :, ", k, ") = ", &
230 pctsrf_t(i, :, k)
231 print *, "phy_ice(", i, ") = ", phy_ice(i)
232 ENDIF
233 END DO
234 ENDDO
235
236 PRINT *, 'Traitement de la sst'
237 call NF95_OPEN('amipbc_sst_1x1.nc', NF90_NOWRITE, ncid)
238
239 call nf95_inq_varid(ncid, "longitude", varid)
240 call nf95_gw_var(ncid, varid, dlon_ini)
241 imdep = size(dlon_ini)
242
243 call nf95_inq_varid(ncid, "latitude", varid)
244 call nf95_gw_var(ncid, varid, dlat_ini)
245 jmdep = size(dlat_ini)
246
247 call nf95_inq_dimid(ncid, "time", dimid)
248 call NF95_INQuire_DIMension(ncid, dimid, nclen=lmdep)
249 print *, 'lmdep = ', lmdep
250 !PM28/02/2002 : nouvelle coord temporelle fichiers AMIP pas en jours
251 ! Ici on suppose qu'on a 12 mois (de 30 jours).
252 IF (lmdep /= 12) stop 'Unknown AMIP file: not 12 months?'
253
254 ALLOCATE(champ(imdep, jmdep), champtime(iim, jjm + 1, lmdep))
255 IF(extrap) THEN
256 ALLOCATE(work(imdep, jmdep))
257 ENDIF
258 ALLOCATE(dlon(imdep), dlat(jmdep))
259 call NF95_INQ_VARID(ncid, 'tosbcs', varid)
260
261 DO l = 1, lmdep
262 call NF95_GET_VAR(ncid, varid, champ, start=(/1, 1, l/))
263 CALL conf_dat2d(dlon_ini, dlat_ini, dlon, dlat, champ)
264 IF (extrap) THEN
265 CALL extrapol(champ, imdep, jmdep, 999999., .TRUE., .TRUE., 2, work)
266 ENDIF
267
268 CALL inter_barxy(dlon, dlat(:jmdep -1), champ, rlonu(:iim), rlatv, &
269 champtime(:, :, l))
270 ENDDO
271
272 call NF95_CLOSE(ncid)
273
274 DEALLOCATE(dlon, dlat, champ, dlon_ini, dlat_ini)
275 allocate(yder(lmdep))
276
277 ! interpolation temporelle
278 DO j = 1, jjm + 1
279 DO i = 1, iim
280 yder = SPLINE(tmidmonth, champtime(i, j, :))
281 DO k = 1, 360
282 champan(i, j, k) = SPLINT(tmidmonth, champtime(i, j, :), yder, &
283 real(k-1))
284 ENDDO
285 ENDDO
286 ENDDO
287
288 deallocate(champtime, yder)
289 champan(iim + 1, :, :) = champan(1, :, :)
290
291 !IM14/03/2002 : SST amipbc greater then 271.38
292 PRINT *, 'SUB. limit_netcdf.F IM : SST Amipbc >= 271.38 '
293 DO k = 1, 360
294 DO j = 1, jjm + 1
295 DO i = 1, iim
296 champan(i, j, k) = amax1(champan(i, j, k), 271.38)
297 ENDDO
298 champan(iim + 1, j, k) = champan(1, j, k)
299 ENDDO
300 ENDDO
301 forall (k = 1:360) phy_sst(:, k) = pack(champan(:, :, k), dyn_phy)
302
303 PRINT *, "Traitement de l'albedo..."
304 call NF95_OPEN('Albedo.nc', NF90_NOWRITE, ncid)
305
306 call nf95_inq_varid(ncid, "longitude", varid)
307 call nf95_gw_var(ncid, varid, dlon_ini)
308 imdep = size(dlon_ini)
309
310 call nf95_inq_varid(ncid, "latitude", varid)
311 call nf95_gw_var(ncid, varid, dlat_ini)
312 jmdep = size(dlat_ini)
313
314 call nf95_inq_varid(ncid, "temps", varid)
315 call nf95_gw_var(ncid, varid, timeyear)
316 lmdep = size(timeyear)
317
318 ALLOCATE(champ(imdep, jmdep), champtime(iim, jjm + 1, lmdep))
319 ALLOCATE(dlon(imdep), dlat(jmdep))
320 call NF95_INQ_VARID(ncid, 'ALBEDO', varid)
321
322 DO l = 1, lmdep
323 PRINT *, "timeyear(", l, ") =", timeyear(l)
324 call NF95_GET_VAR(ncid, varid, champ, start=(/1, 1, l/))
325 CALL conf_dat2d(dlon_ini, dlat_ini, dlon, dlat, champ)
326 CALL inter_barxy(dlon, dlat(:jmdep-1), champ, rlonu(:iim), rlatv, &
327 champtime(:, :, l))
328 ENDDO
329
330 call NF95_CLOSE(ncid)
331
332 deallocate(dlon_ini, dlat_ini)
333 allocate(yder(lmdep))
334
335 ! interpolation temporelle
336 DO j = 1, jjm + 1
337 DO i = 1, iim
338 yder = SPLINE(timeyear, champtime(i, j, :))
339 DO k = 1, 360
340 champan(i, j, k) = SPLINT(timeyear, champtime(i, j, :), yder, &
341 real(k-1))
342 ENDDO
343 ENDDO
344 ENDDO
345 deallocate(timeyear)
346
347 champan(iim + 1, :, :) = champan(1, :, :)
348 forall (k = 1:360) phy_alb(:, k) = pack(champan(:, :, k), dyn_phy)
349
350 DO k = 1, 360
351 DO i = 1, klon
352 phy_bil(i, k) = 0.0
353 ENDDO
354 ENDDO
355
356 PRINT *, 'Ecriture du fichier limit.nc...'
357
358 call NF95_CREATE("limit.nc", NF90_CLOBBER, nid)
359
360 call NF95_PUT_ATT(nid, NF90_GLOBAL, "title", &
361 "Fichier conditions aux limites")
362 call NF95_DEF_DIM(nid, "points_physiques", klon, ndim)
363 call NF95_DEF_DIM(nid, "time", NF90_UNLIMITED, ntim)
364
365 call NF95_DEF_VAR(nid, "TEMPS", NF90_FLOAT, ntim, id_tim)
366 call NF95_PUT_ATT(nid, id_tim, "title", "Jour dans l annee")
367
368 call NF95_DEF_VAR(nid, "FOCE", NF90_FLOAT, dimids=(/ndim, ntim/), &
369 varid=id_foce)
370 call NF95_PUT_ATT(nid, id_FOCE, "title", "Fraction ocean")
371
372 call NF95_DEF_VAR(nid, "FSIC", NF90_FLOAT, dimids=(/ndim, ntim/), &
373 varid=id_FSIC)
374 call NF95_PUT_ATT(nid, id_FSIC, "title", "Fraction glace de mer")
375
376 call NF95_DEF_VAR(nid, "FTER", NF90_FLOAT, dimids=(/ndim, ntim/), &
377 varid=id_FTER)
378 call NF95_PUT_ATT(nid, id_FTER, "title", "Fraction terre")
379
380 call NF95_DEF_VAR(nid, "FLIC", NF90_FLOAT, dimids=(/ndim, ntim/), &
381 varid=id_FLIC)
382 call NF95_PUT_ATT(nid, id_FLIC, "title", "Fraction land ice")
383
384 call NF95_DEF_VAR(nid, "SST", NF90_FLOAT, dimids=(/ndim, ntim/), &
385 varid=id_SST)
386 call NF95_PUT_ATT(nid, id_SST, "title", &
387 "Temperature superficielle de la mer")
388
389 call NF95_DEF_VAR(nid, "BILS", NF90_FLOAT, dimids=(/ndim, ntim/), &
390 varid=id_BILS)
391 call NF95_PUT_ATT(nid, id_BILS, "title", "Reference flux de chaleur au sol")
392
393 call NF95_DEF_VAR(nid, "ALB", NF90_FLOAT, dimids=(/ndim, ntim/), &
394 varid=id_ALB)
395 call NF95_PUT_ATT(nid, id_ALB, "title", "Albedo a la surface")
396
397 call NF95_DEF_VAR(nid, "RUG", NF90_FLOAT, dimids=(/ndim, ntim/), &
398 varid=id_RUG)
399 call NF95_PUT_ATT(nid, id_RUG, "title", "Rugosite")
400
401 call NF95_ENDDEF(nid)
402
403 DO k = 1, 360
404 call NF95_PUT_VAR(nid, id_tim, REAL(k), (/k/))
405 call NF95_PUT_VAR(nid, id_FOCE, pctsrf_t(:, is_oce, k), start=(/1, k/))
406 call NF95_PUT_VAR(nid, id_FSIC, pctsrf_t(:, is_sic, k), start=(/1, k/))
407 call NF95_PUT_VAR(nid, id_FTER, pctsrf_t(:, is_ter, k), start=(/1, k/))
408 call NF95_PUT_VAR(nid, id_FLIC, pctsrf_t(:, is_lic, k), start=(/1, k/))
409 call NF95_PUT_VAR(nid, id_SST, phy_sst(:, k), start=(/1, k/))
410 call NF95_PUT_VAR(nid, id_BILS, phy_bil(:, k), start=(/1, k/))
411 call NF95_PUT_VAR(nid, id_ALB, phy_alb(:, k), start=(/1, k/))
412 call NF95_PUT_VAR(nid, id_RUG, phy_rug(:, k), start=(/1, k/))
413 ENDDO
414
415 call NF95_CLOSE(nid)
416
417 END SUBROUTINE limit
418
419 end module limit_mod

  ViewVC Help
Powered by ViewVC 1.1.21