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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 57 - (show annotations)
Mon Jan 30 12:54:02 2012 UTC (12 years, 3 months ago) by guez
File size: 14184 byte(s)
Write used namelists to file "" instead of standard output.

Avoid aliasing in "inidissip" in calls to "divgrad2", "divgrad",
"gradiv2", "gradiv", "nxgraro2" and "nxgrarot". Add a degenerate
dimension to arrays so they have rank 3, like the dummy arguments in
"divgrad2", "divgrad", "gradiv2", "gradiv", "nxgraro2" and "nxgrarot".

Extract the initialization part from "bilan_dyn" and make a separate
procedure, "init_dynzon", from it.

Move variables from modules "iniprint" and "logic" to module
"conf_gcm_m".

Promote internal procedures of "fxy" to private procedures of module
"fxy_m".

Extracted documentation from "inigeom". Removed useless "save"
attributes. Removed useless intermediate variables. Extracted
processing of poles from loop on latitudes. Write coordinates to file
"longitude_latitude.txt" instead of standard output.

Do not use ozone tracer for radiative transfer.

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

  ViewVC Help
Powered by ViewVC 1.1.21