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

Annotation of /trunk/Sources/dyn3d/limit.f

Parent Directory Parent Directory | Revision Log Revision Log


Revision 139 - (hide annotations)
Tue May 26 17:46:03 2015 UTC (8 years, 11 months ago) by guez
File size: 13921 byte(s)
dynetat0 read rlonu, rlatu, rlonv, rlatv, cu_2d, cv_2d, aire_2d from
"start.nc" and then these variables were overwritten by
inigeom. Corrected this. Now, inigeom does not compute rlonu, rlatu,
rlonv and rlatv. Moreover, cu_2d, cv_2d, aire_2d are not written to
"restart.nc". Since xprimu, xprimv, xprimm025, xprimp025, rlatu1,
rlatu2, yprimu1, yprimu2 are computed at the same time as rlonu,
rlatu, rlonv, rlatv, and since it would not be convenient to separate
those computations, we decide to write xprimu, xprimv, xprimm025,
xprimp025, rlatu1, rlatu2, yprimu1, yprimu2 into "restart.nc", read
them from "start.nc" and not compute them in inigeom. So, in summary,
"start.nc" contains all the coordinates and their derivatives, and
inigeom only computes the 2D-variables.

Technical details:

Moved variables rlatu, rlonv, rlonu, rlatv, xprimu, xprimv from module
comgeom to module dynetat0_m. Upgraded local variables rlatu1,
yprimu1, rlatu2, yprimu2, xprimm025, xprimp025 of procedure inigeom to
variables of module dynetat0_m.

Removed unused local variable yprimu of procedure inigeom and
corresponding argument yyprimu of fyhyp.

Moved variables clat, clon, grossismx, grossismy, dzoomx, dzoomy,
taux, tauy from module serre to module dynetat0_m (since they are read
from "start.nc"). The default values are now defined in read_serre
instead of in the declarations. Changed name of module serre to
read_serre_m, no more module variable here.

The calls to fxhyp and fyhyp are moved from inigeom to etat0.

Side effects in programs other than gcm: etat0 and read_serre write
variables of module dynetat0; the programs test_fxyp and
test_inter_barxy need more source files.

Removed unused arguments len and nd of cv3_tracer. Removed unused
argument PPSOL of LWU.

Bug fix in test_inter_barxy: forgotten call to read_serre.

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 conf_dat2d_m, only: conf_dat2d
16 guez 3 use dimens_m, only: iim, jjm
17     use dimphy, only: klon, zmasq
18 guez 139 use dynetat0_m, only: rlonu, rlatv
19 guez 3 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 98 use netcdf95, only: NF95_CLOSE, NF95_CREATE, NF95_DEF_DIM, nf95_def_var, &
24 guez 97 nf95_enddef, nf95_get_var, nf95_gw_var, nf95_inq_dimid, &
25 guez 98 nf95_inq_varid, nf95_inquire_dimension, NF95_OPEN, NF95_PUT_ATT, &
26     NF95_PUT_VAR
27     use netcdf, only: NF90_CLOBBER, NF90_FLOAT, NF90_GLOBAL, NF90_NOWRITE, &
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 guez 139 ! (extrapolation de donn\'ees, comme pour les SST lorsque le fichier
37     ! ne contient pas uniquement des points oc\'eaniques)
38 guez 3
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 guez 139 ! Pour le champ de d\'epart:
48 guez 3 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 guez 139 ! Pour le champ interpol\'e 3D :
56 guez 3 REAL, allocatable:: champtime(:, :, :)
57     REAL champan(iim + 1, jjm + 1, 360)
58    
59     ! Pour l'inteprolation verticale :
60     REAL, allocatable:: yder(:)
61    
62     INTEGER nid, ndim, ntim
63     INTEGER id_tim
64     INTEGER id_SST, id_BILS, id_RUG, id_ALB
65     INTEGER id_FOCE, id_FSIC, id_FTER, id_FLIC
66    
67     INTEGER i, j, k, l
68     INTEGER ncid, varid, dimid
69    
70     REAL, parameter:: tmidmonth(12) = (/(15. + 30. * i, i = 0, 11)/)
71    
72     namelist /limit_nml/extrap
73    
74     !--------------------
75    
76     print *, "Call sequence information: limit"
77    
78     print *, "Enter namelist 'limit_nml'."
79 guez 98 read(unit=*, nml=limit_nml)
80 guez 57 write(unit_nml, nml=limit_nml)
81 guez 3
82 guez 68 PRINT *, 'Processing rugosity...'
83 guez 3
84     call NF95_OPEN('Rugos.nc', NF90_NOWRITE, ncid)
85    
86 guez 15 ! Read coordinate variables:
87    
88 guez 22 call nf95_inq_varid(ncid, "longitude", varid)
89     call nf95_gw_var(ncid, varid, dlon_ini)
90 guez 3 imdep = size(dlon_ini)
91    
92 guez 22 call nf95_inq_varid(ncid, "latitude", varid)
93     call nf95_gw_var(ncid, varid, dlat_ini)
94 guez 3 jmdep = size(dlat_ini)
95    
96 guez 22 call nf95_inq_varid(ncid, "temps", varid)
97     call nf95_gw_var(ncid, varid, timeyear)
98 guez 3 lmdep = size(timeyear)
99    
100     ALLOCATE(champ(imdep, jmdep), champtime(iim, jjm + 1, lmdep))
101     allocate(dlon(imdep), dlat(jmdep))
102     call NF95_INQ_VARID(ncid, 'RUGOS', varid)
103    
104 guez 15 ! Read the primary variable day by day and regrid horizontally,
105     ! result in "champtime":
106 guez 3 DO l = 1, lmdep
107 guez 97 call NF95_GET_VAR(ncid, varid, champ, start=(/1, 1, l/))
108 guez 3 CALL conf_dat2d(dlon_ini, dlat_ini, dlon, dlat, champ)
109     CALL inter_barxy(dlon, dlat(:jmdep -1), LOG(champ), rlonu(:iim), &
110     rlatv, champtime(:, :, l))
111     champtime(:, :, l) = EXP(champtime(:, :, l))
112 guez 15 where (nint(mask(:iim, :)) /= 1) champtime(:, :, l) = 0.001
113 guez 3 end do
114    
115     call NF95_CLOSE(ncid)
116    
117     DEALLOCATE(dlon, dlat, champ, dlon_ini, dlat_ini)
118     allocate(yder(lmdep))
119    
120 guez 15 ! Interpolate monthly values to daily values, at each horizontal position:
121 guez 3 DO j = 1, jjm + 1
122     DO i = 1, iim
123 guez 68 yder = SPLINE(timeyear, champtime(i, j, :))
124 guez 3 DO k = 1, 360
125     champan(i, j, k) = SPLINT(timeyear, champtime(i, j, :), yder, &
126     real(k-1))
127     ENDDO
128     ENDDO
129     ENDDO
130    
131     deallocate(timeyear, champtime, yder)
132     champan(iim + 1, :, :) = champan(1, :, :)
133     forall (k = 1:360) phy_rug(:, k) = pack(champan(:, :, k), dyn_phy)
134    
135     ! Process sea ice:
136    
137     PRINT *, 'Processing sea ice...'
138     call NF95_OPEN('amipbc_sic_1x1.nc', NF90_NOWRITE, ncid)
139    
140 guez 22 call nf95_inq_varid(ncid, "longitude", varid)
141     call nf95_gw_var(ncid, varid, dlon_ini)
142 guez 3 imdep = size(dlon_ini)
143    
144 guez 22 call nf95_inq_varid(ncid, "latitude", varid)
145     call nf95_gw_var(ncid, varid, dlat_ini)
146 guez 3 jmdep = size(dlat_ini)
147    
148     call nf95_inq_dimid(ncid, "time", dimid)
149 guez 34 call NF95_INQuire_DIMension(ncid, dimid, nclen=lmdep)
150 guez 3 print *, 'lmdep = ', lmdep
151 guez 139 ! Coordonn\'ee temporelle fichiers AMIP pas en jours. Ici on suppose
152 guez 68 ! qu'on a 12 mois (de 30 jours).
153     IF (lmdep /= 12) then
154     print *, 'Unknown AMIP file: not 12 months?'
155     STOP 1
156     end IF
157 guez 3
158     ALLOCATE(champ(imdep, jmdep), champtime(iim, jjm + 1, lmdep))
159 guez 98 ALLOCATE(dlon(imdep), dlat(jmdep))
160 guez 3 call NF95_INQ_VARID(ncid, 'sicbcs', varid)
161     DO l = 1, lmdep
162 guez 97 call NF95_GET_VAR(ncid, varid, champ, start=(/1, 1, l/))
163 guez 3 CALL conf_dat2d(dlon_ini, dlat_ini, dlon, dlat, champ)
164 guez 98 CALL inter_barxy(dlon, dlat(:jmdep -1), champ, rlonu(:iim), rlatv, &
165 guez 3 champtime(:, :, l))
166     ENDDO
167    
168     call NF95_CLOSE(ncid)
169    
170     DEALLOCATE(dlon, dlat, champ, dlon_ini, dlat_ini)
171     PRINT *, 'Interpolation temporelle'
172     allocate(yder(lmdep))
173    
174     DO j = 1, jjm + 1
175     DO i = 1, iim
176 guez 68 yder = SPLINE(tmidmonth, champtime(i, j, :))
177 guez 3 DO k = 1, 360
178     champan(i, j, k) = SPLINT(tmidmonth, champtime(i, j, :), yder, &
179     real(k-1))
180     ENDDO
181     ENDDO
182     ENDDO
183    
184     deallocate(champtime, yder)
185    
186     ! Convert from percentage to normal fraction and keep sea ice
187     ! between 0 and 1:
188     champan(:iim, :, :) = max(0., (min(1., (champan(:iim, :, :) / 100.))))
189     champan(iim + 1, :, :) = champan(1, :, :)
190    
191     DO k = 1, 360
192 guez 68 phy_ice = pack(champan(:, :, k), dyn_phy)
193 guez 3
194     ! (utilisation de la sous-maille fractionnelle tandis que l'ancien
195     ! codage utilisait l'indicateur du sol (0, 1, 2, 3))
196     ! PB en attendant de mettre fraction de terre
197 guez 68 WHERE (phy_ice < EPSFRA) phy_ice = 0.
198 guez 3
199     pctsrf_t(:, is_ter, k) = pctsrf(:, is_ter)
200     pctsrf_t(:, is_lic, k) = pctsrf(:, is_lic)
201 guez 68 pctsrf_t(:, is_sic, k) = max(phy_ice - pctsrf_t(:, is_lic, k), 0.)
202 guez 139 ! Il y a des cas o\`u il y a de la glace dans landiceref et
203 guez 3 ! pas dans AMIP
204 guez 68 WHERE (1. - zmasq < EPSFRA)
205 guez 3 pctsrf_t(:, is_sic, k) = 0.
206     pctsrf_t(:, is_oce, k) = 0.
207     elsewhere
208 guez 68 where (pctsrf_t(:, is_sic, k) >= 1 - zmasq)
209     pctsrf_t(:, is_sic, k) = 1. - zmasq
210 guez 3 pctsrf_t(:, is_oce, k) = 0.
211     ELSEwhere
212 guez 68 pctsrf_t(:, is_oce, k) = 1. - zmasq - pctsrf_t(:, is_sic, k)
213 guez 3 where (pctsrf_t(:, is_oce, k) < EPSFRA)
214     pctsrf_t(:, is_oce, k) = 0.
215 guez 68 pctsrf_t(:, is_sic, k) = 1 - zmasq
216 guez 3 end where
217     end where
218     end where
219    
220     DO i = 1, klon
221     if (pctsrf_t(i, is_oce, k) < 0.) then
222 guez 139 print *, 'Bad surface fraction: pctsrf_t(', i, &
223 guez 3 ', is_oce, ', k, ') = ', pctsrf_t(i, is_oce, k)
224     ENDIF
225     IF (abs(pctsrf_t(i, is_ter, k) + pctsrf_t(i, is_lic, k) &
226     + pctsrf_t(i, is_oce, k) + pctsrf_t(i, is_sic, k) - 1.) &
227     > EPSFRA) THEN
228 guez 139 print *, 'Bad surface fraction:'
229 guez 3 print *, "pctsrf_t(", i, ", :, ", k, ") = ", &
230     pctsrf_t(i, :, k)
231     print *, "phy_ice(", i, ") = ", phy_ice(i)
232     ENDIF
233     END DO
234     ENDDO
235    
236     PRINT *, 'Traitement de la sst'
237     call NF95_OPEN('amipbc_sst_1x1.nc', NF90_NOWRITE, ncid)
238    
239 guez 22 call nf95_inq_varid(ncid, "longitude", varid)
240     call nf95_gw_var(ncid, varid, dlon_ini)
241 guez 3 imdep = size(dlon_ini)
242    
243 guez 22 call nf95_inq_varid(ncid, "latitude", varid)
244     call nf95_gw_var(ncid, varid, dlat_ini)
245 guez 3 jmdep = size(dlat_ini)
246    
247     call nf95_inq_dimid(ncid, "time", dimid)
248 guez 34 call NF95_INQuire_DIMension(ncid, dimid, nclen=lmdep)
249 guez 3 print *, 'lmdep = ', lmdep
250     !PM28/02/2002 : nouvelle coord temporelle fichiers AMIP pas en jours
251     ! Ici on suppose qu'on a 12 mois (de 30 jours).
252     IF (lmdep /= 12) stop 'Unknown AMIP file: not 12 months?'
253    
254 guez 68 ALLOCATE(champ(imdep, jmdep), champtime(iim, jjm + 1, lmdep))
255     IF(extrap) THEN
256 guez 98 ALLOCATE(work(imdep, jmdep))
257 guez 3 ENDIF
258 guez 68 ALLOCATE(dlon(imdep), dlat(jmdep))
259 guez 3 call NF95_INQ_VARID(ncid, 'tosbcs', varid)
260    
261     DO l = 1, lmdep
262 guez 97 call NF95_GET_VAR(ncid, varid, champ, start=(/1, 1, l/))
263 guez 3 CALL conf_dat2d(dlon_ini, dlat_ini, dlon, dlat, champ)
264 guez 68 IF (extrap) THEN
265 guez 3 CALL extrapol(champ, imdep, jmdep, 999999., .TRUE., .TRUE., 2, work)
266     ENDIF
267    
268 guez 98 CALL inter_barxy(dlon, dlat(:jmdep -1), champ, rlonu(:iim), rlatv, &
269 guez 68 champtime(:, :, l))
270 guez 3 ENDDO
271    
272     call NF95_CLOSE(ncid)
273    
274     DEALLOCATE(dlon, dlat, champ, dlon_ini, dlat_ini)
275     allocate(yder(lmdep))
276    
277     ! interpolation temporelle
278     DO j = 1, jjm + 1
279     DO i = 1, iim
280 guez 68 yder = SPLINE(tmidmonth, champtime(i, j, :))
281 guez 3 DO k = 1, 360
282     champan(i, j, k) = SPLINT(tmidmonth, champtime(i, j, :), yder, &
283     real(k-1))
284     ENDDO
285     ENDDO
286     ENDDO
287    
288     deallocate(champtime, yder)
289     champan(iim + 1, :, :) = champan(1, :, :)
290    
291     !IM14/03/2002 : SST amipbc greater then 271.38
292     PRINT *, 'SUB. limit_netcdf.F IM : SST Amipbc >= 271.38 '
293     DO k = 1, 360
294     DO j = 1, jjm + 1
295     DO i = 1, iim
296     champan(i, j, k) = amax1(champan(i, j, k), 271.38)
297     ENDDO
298     champan(iim + 1, j, k) = champan(1, j, k)
299     ENDDO
300     ENDDO
301     forall (k = 1:360) phy_sst(:, k) = pack(champan(:, :, k), dyn_phy)
302    
303 guez 97 PRINT *, "Traitement de l'albedo..."
304 guez 3 call NF95_OPEN('Albedo.nc', NF90_NOWRITE, ncid)
305    
306 guez 22 call nf95_inq_varid(ncid, "longitude", varid)
307     call nf95_gw_var(ncid, varid, dlon_ini)
308 guez 3 imdep = size(dlon_ini)
309    
310 guez 22 call nf95_inq_varid(ncid, "latitude", varid)
311     call nf95_gw_var(ncid, varid, dlat_ini)
312 guez 3 jmdep = size(dlat_ini)
313    
314 guez 22 call nf95_inq_varid(ncid, "temps", varid)
315     call nf95_gw_var(ncid, varid, timeyear)
316 guez 3 lmdep = size(timeyear)
317    
318 guez 98 ALLOCATE(champ(imdep, jmdep), champtime(iim, jjm + 1, lmdep))
319     ALLOCATE(dlon(imdep), dlat(jmdep))
320 guez 3 call NF95_INQ_VARID(ncid, 'ALBEDO', varid)
321    
322     DO l = 1, lmdep
323 guez 97 PRINT *, "timeyear(", l, ") =", timeyear(l)
324     call NF95_GET_VAR(ncid, varid, champ, start=(/1, 1, l/))
325 guez 3 CALL conf_dat2d(dlon_ini, dlat_ini, dlon, dlat, champ)
326     CALL inter_barxy(dlon, dlat(:jmdep-1), champ, rlonu(:iim), rlatv, &
327 guez 68 champtime(:, :, l))
328 guez 3 ENDDO
329    
330     call NF95_CLOSE(ncid)
331    
332     deallocate(dlon_ini, dlat_ini)
333     allocate(yder(lmdep))
334    
335     ! interpolation temporelle
336     DO j = 1, jjm + 1
337     DO i = 1, iim
338 guez 68 yder = SPLINE(timeyear, champtime(i, j, :))
339 guez 3 DO k = 1, 360
340     champan(i, j, k) = SPLINT(timeyear, champtime(i, j, :), yder, &
341     real(k-1))
342     ENDDO
343     ENDDO
344     ENDDO
345     deallocate(timeyear)
346    
347     champan(iim + 1, :, :) = champan(1, :, :)
348     forall (k = 1:360) phy_alb(:, k) = pack(champan(:, :, k), dyn_phy)
349    
350     DO k = 1, 360
351     DO i = 1, klon
352     phy_bil(i, k) = 0.0
353     ENDDO
354     ENDDO
355    
356 guez 98 PRINT *, 'Ecriture du fichier limit.nc...'
357 guez 3
358     call NF95_CREATE("limit.nc", NF90_CLOBBER, nid)
359    
360 guez 98 call NF95_PUT_ATT(nid, NF90_GLOBAL, "title", &
361 guez 3 "Fichier conditions aux limites")
362 guez 98 call NF95_DEF_DIM(nid, "points_physiques", klon, ndim)
363     call NF95_DEF_DIM(nid, "time", NF90_UNLIMITED, ntim)
364 guez 3
365 guez 98 call NF95_DEF_VAR(nid, "TEMPS", NF90_FLOAT, ntim, id_tim)
366     call NF95_PUT_ATT(nid, id_tim, "title", "Jour dans l annee")
367 guez 3
368 guez 98 call NF95_DEF_VAR(nid, "FOCE", NF90_FLOAT, dimids=(/ndim, ntim/), &
369     varid=id_foce)
370     call NF95_PUT_ATT(nid, id_FOCE, "title", "Fraction ocean")
371 guez 3
372 guez 98 call NF95_DEF_VAR(nid, "FSIC", NF90_FLOAT, dimids=(/ndim, ntim/), &
373     varid=id_FSIC)
374     call NF95_PUT_ATT(nid, id_FSIC, "title", "Fraction glace de mer")
375 guez 3
376 guez 98 call NF95_DEF_VAR(nid, "FTER", NF90_FLOAT, dimids=(/ndim, ntim/), &
377     varid=id_FTER)
378     call NF95_PUT_ATT(nid, id_FTER, "title", "Fraction terre")
379 guez 3
380 guez 98 call NF95_DEF_VAR(nid, "FLIC", NF90_FLOAT, dimids=(/ndim, ntim/), &
381     varid=id_FLIC)
382     call NF95_PUT_ATT(nid, id_FLIC, "title", "Fraction land ice")
383 guez 3
384 guez 98 call NF95_DEF_VAR(nid, "SST", NF90_FLOAT, dimids=(/ndim, ntim/), &
385     varid=id_SST)
386     call NF95_PUT_ATT(nid, id_SST, "title", &
387 guez 3 "Temperature superficielle de la mer")
388    
389 guez 98 call NF95_DEF_VAR(nid, "BILS", NF90_FLOAT, dimids=(/ndim, ntim/), &
390     varid=id_BILS)
391     call NF95_PUT_ATT(nid, id_BILS, "title", "Reference flux de chaleur au sol")
392    
393     call NF95_DEF_VAR(nid, "ALB", NF90_FLOAT, dimids=(/ndim, ntim/), &
394     varid=id_ALB)
395     call NF95_PUT_ATT(nid, id_ALB, "title", "Albedo a la surface")
396    
397     call NF95_DEF_VAR(nid, "RUG", NF90_FLOAT, dimids=(/ndim, ntim/), &
398     varid=id_RUG)
399     call NF95_PUT_ATT(nid, id_RUG, "title", "Rugosite")
400    
401 guez 3 call NF95_ENDDEF(nid)
402    
403     DO k = 1, 360
404 guez 98 call NF95_PUT_VAR(nid, id_tim, REAL(k), (/k/))
405     call NF95_PUT_VAR(nid, id_FOCE, pctsrf_t(:, is_oce, k), start=(/1, k/))
406     call NF95_PUT_VAR(nid, id_FSIC, pctsrf_t(:, is_sic, k), start=(/1, k/))
407     call NF95_PUT_VAR(nid, id_FTER, pctsrf_t(:, is_ter, k), start=(/1, k/))
408     call NF95_PUT_VAR(nid, id_FLIC, pctsrf_t(:, is_lic, k), start=(/1, k/))
409     call NF95_PUT_VAR(nid, id_SST, phy_sst(:, k), start=(/1, k/))
410     call NF95_PUT_VAR(nid, id_BILS, phy_bil(:, k), start=(/1, k/))
411     call NF95_PUT_VAR(nid, id_ALB, phy_alb(:, k), start=(/1, k/))
412     call NF95_PUT_VAR(nid, id_RUG, phy_rug(:, k), start=(/1, k/))
413 guez 3 ENDDO
414    
415     call NF95_CLOSE(nid)
416    
417     END SUBROUTINE limit
418    
419     end module limit_mod

  ViewVC Help
Powered by ViewVC 1.1.21