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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 57 - (hide annotations)
Mon Jan 30 12:54:02 2012 UTC (12 years, 4 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 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 57 use numer_rec, only: spline, splint
30     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     ! 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 guez 34 call NF95_INQuire_DIMension(ncid, dimid, nclen=lmdep)
156 guez 3 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 guez 34 call NF95_INQuire_DIMension(ncid, dimid, nclen=lmdep)
255 guez 3 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