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

Annotation of /trunk/dyn3d/limit.f

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.21