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

Annotation of /trunk/dyn3d/limit.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 25 - (hide annotations)
Fri Mar 5 16:43:45 2010 UTC (14 years, 2 months ago) by guez
Original Path: trunk/libf/dyn3d/limit.f90
File size: 14144 byte(s)
Simplified "etat0_lim.sh" and "gcm.sh" because the full versions
depended on personal arrangements for directories and machines.

Translated included files into modules. Encapsulated procedures into modules.

Moved variables from module "comgeom" to local variables of
"inigeom". Deleted some unused variables in "comgeom".

Moved variable "day_ini" from module "temps" to module "dynetat0_m".

Removed useless test on variable "time" and useless "close" statement
in procedure "leapfrog".

Removed useless call to "inigeom" in procedure "limit".

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     use dimens_m, only: iim, jjm
16     use indicesol, only: epsfra, nbsrf, is_ter, is_oce, is_lic, is_sic
17     use dimphy, only: klon, zmasq
18     use comgeom, only: rlonu, rlatv
19     use etat0_mod, only: pctsrf
20 guez 15 use start_init_orog_m, only: mask
21 guez 3 use conf_dat2d_m, only: conf_dat2d
22     use inter_barxy_m, only: inter_barxy
23 guez 13 use numer_rec, only: spline, splint
24 guez 3 use grid_change, only: dyn_phy
25    
26 guez 22 use netcdf95, only: handle_err, nf95_gw_var, NF95_CLOSE, NF95_DEF_DIM, &
27 guez 7 nf95_enddef, NF95_CREATE, nf95_inq_dimid, nf95_inquire_dimension, &
28     nf95_inq_varid, NF95_OPEN
29     use netcdf, only: NF90_CLOBBER, nf90_def_var, NF90_FLOAT, NF90_GET_VAR, &
30     NF90_GLOBAL, NF90_NOWRITE, NF90_PUT_ATT, NF90_PUT_VAR, &
31     NF90_UNLIMITED
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     write(unit=*, nml=limit_nml)
84    
85     ! Process rugosity:
86    
87     PRINT *, 'Processing rugosity...'
88     call NF95_OPEN('Rugos.nc', NF90_NOWRITE, ncid)
89    
90 guez 15 ! Read coordinate variables:
91    
92 guez 22 call nf95_inq_varid(ncid, "longitude", varid)
93     call nf95_gw_var(ncid, varid, dlon_ini)
94 guez 3 imdep = size(dlon_ini)
95    
96 guez 22 call nf95_inq_varid(ncid, "latitude", varid)
97     call nf95_gw_var(ncid, varid, dlat_ini)
98 guez 3 jmdep = size(dlat_ini)
99    
100 guez 22 call nf95_inq_varid(ncid, "temps", varid)
101     call nf95_gw_var(ncid, varid, timeyear)
102 guez 3 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 guez 15 ! Read the primary variable day by day and regrid horizontally,
109     ! result in "champtime":
110 guez 3 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 guez 15 where (nint(mask(:iim, :)) /= 1) champtime(:, :, l) = 0.001
119 guez 3 end do
120    
121     call NF95_CLOSE(ncid)
122    
123     DEALLOCATE(dlon, dlat, champ, dlon_ini, dlat_ini)
124     allocate(yder(lmdep))
125    
126 guez 15 ! Interpolate monthly values to daily values, at each horizontal position:
127 guez 3 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 guez 22 call nf95_inq_varid(ncid, "longitude", varid)
147     call nf95_gw_var(ncid, varid, dlon_ini)
148 guez 3 imdep = size(dlon_ini)
149    
150 guez 22 call nf95_inq_varid(ncid, "latitude", varid)
151     call nf95_gw_var(ncid, varid, dlat_ini)
152 guez 3 jmdep = size(dlat_ini)
153    
154     call nf95_inq_dimid(ncid, "time", dimid)
155     call NF95_INQuire_DIMension(ncid, dimid, len=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 guez 22 call nf95_inq_varid(ncid, "longitude", varid)
246     call nf95_gw_var(ncid, varid, dlon_ini)
247 guez 3 imdep = size(dlon_ini)
248    
249 guez 22 call nf95_inq_varid(ncid, "latitude", varid)
250     call nf95_gw_var(ncid, varid, dlat_ini)
251 guez 3 jmdep = size(dlat_ini)
252    
253     call nf95_inq_dimid(ncid, "time", dimid)
254     call NF95_INQuire_DIMension(ncid, dimid, len=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 guez 22 call nf95_inq_varid(ncid, "longitude", varid)
317     call nf95_gw_var(ncid, varid, dlon_ini)
318 guez 3 imdep = size(dlon_ini)
319    
320 guez 22 call nf95_inq_varid(ncid, "latitude", varid)
321     call nf95_gw_var(ncid, varid, dlat_ini)
322 guez 3 jmdep = size(dlat_ini)
323    
324 guez 22 call nf95_inq_varid(ncid, "temps", varid)
325     call nf95_gw_var(ncid, varid, timeyear)
326 guez 3 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