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

Contents of /trunk/libf/dyn3d/limit.f90

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.21