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

Annotation of /trunk/dyn3d/limit.f

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.21