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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 61 - (show annotations)
Fri Apr 20 14:58:43 2012 UTC (12 years ago) by guez
File size: 14187 byte(s)
No more included file in LMDZE, not even "netcdf.inc".

Created a variable containing the list of common source files in
GNUmakefile. So we now also see clearly files that are specific to
each program.

Split module "histcom". Assembled resulting files in directory
"Histcom".

Removed aliasing in calls to "laplacien".

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_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 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