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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 139 - (show 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 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 conf_dat2d_m, only: conf_dat2d
16 use dimens_m, only: iim, jjm
17 use dimphy, only: klon, zmasq
18 use dynetat0_m, only: rlonu, rlatv
19 use etat0_mod, only: pctsrf
20 use grid_change, only: dyn_phy
21 use indicesol, only: epsfra, nbsrf, is_ter, is_oce, is_lic, is_sic
22 use inter_barxy_m, only: inter_barxy
23 use netcdf95, only: NF95_CLOSE, NF95_CREATE, NF95_DEF_DIM, nf95_def_var, &
24 nf95_enddef, nf95_get_var, nf95_gw_var, nf95_inq_dimid, &
25 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 use numer_rec_95, only: spline, splint
30 use start_init_orog_m, only: mask
31 use unit_nml_m, only: unit_nml
32
33 ! Variables local to the procedure:
34
35 LOGICAL:: extrap = .FALSE.
36 ! (extrapolation de donn\'ees, comme pour les SST lorsque le fichier
37 ! ne contient pas uniquement des points oc\'eaniques)
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\'epart:
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\'e 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 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 read(unit=*, nml=limit_nml)
80 write(unit_nml, nml=limit_nml)
81
82 PRINT *, 'Processing rugosity...'
83
84 call NF95_OPEN('Rugos.nc', NF90_NOWRITE, ncid)
85
86 ! Read coordinate variables:
87
88 call nf95_inq_varid(ncid, "longitude", varid)
89 call nf95_gw_var(ncid, varid, dlon_ini)
90 imdep = size(dlon_ini)
91
92 call nf95_inq_varid(ncid, "latitude", varid)
93 call nf95_gw_var(ncid, varid, dlat_ini)
94 jmdep = size(dlat_ini)
95
96 call nf95_inq_varid(ncid, "temps", varid)
97 call nf95_gw_var(ncid, varid, timeyear)
98 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 ! Read the primary variable day by day and regrid horizontally,
105 ! result in "champtime":
106 DO l = 1, lmdep
107 call NF95_GET_VAR(ncid, varid, champ, start=(/1, 1, l/))
108 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 where (nint(mask(:iim, :)) /= 1) champtime(:, :, l) = 0.001
113 end do
114
115 call NF95_CLOSE(ncid)
116
117 DEALLOCATE(dlon, dlat, champ, dlon_ini, dlat_ini)
118 allocate(yder(lmdep))
119
120 ! Interpolate monthly values to daily values, at each horizontal position:
121 DO j = 1, jjm + 1
122 DO i = 1, iim
123 yder = SPLINE(timeyear, champtime(i, j, :))
124 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 call nf95_inq_varid(ncid, "longitude", varid)
141 call nf95_gw_var(ncid, varid, dlon_ini)
142 imdep = size(dlon_ini)
143
144 call nf95_inq_varid(ncid, "latitude", varid)
145 call nf95_gw_var(ncid, varid, dlat_ini)
146 jmdep = size(dlat_ini)
147
148 call nf95_inq_dimid(ncid, "time", dimid)
149 call NF95_INQuire_DIMension(ncid, dimid, nclen=lmdep)
150 print *, 'lmdep = ', lmdep
151 ! Coordonn\'ee temporelle fichiers AMIP pas en jours. Ici on suppose
152 ! 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
158 ALLOCATE(champ(imdep, jmdep), champtime(iim, jjm + 1, lmdep))
159 ALLOCATE(dlon(imdep), dlat(jmdep))
160 call NF95_INQ_VARID(ncid, 'sicbcs', varid)
161 DO l = 1, lmdep
162 call NF95_GET_VAR(ncid, varid, champ, start=(/1, 1, l/))
163 CALL conf_dat2d(dlon_ini, dlat_ini, dlon, dlat, champ)
164 CALL inter_barxy(dlon, dlat(:jmdep -1), champ, rlonu(:iim), rlatv, &
165 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 yder = SPLINE(tmidmonth, champtime(i, j, :))
177 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 phy_ice = pack(champan(:, :, k), dyn_phy)
193
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 WHERE (phy_ice < EPSFRA) phy_ice = 0.
198
199 pctsrf_t(:, is_ter, k) = pctsrf(:, is_ter)
200 pctsrf_t(:, is_lic, k) = pctsrf(:, is_lic)
201 pctsrf_t(:, is_sic, k) = max(phy_ice - pctsrf_t(:, is_lic, k), 0.)
202 ! Il y a des cas o\`u il y a de la glace dans landiceref et
203 ! pas dans AMIP
204 WHERE (1. - zmasq < EPSFRA)
205 pctsrf_t(:, is_sic, k) = 0.
206 pctsrf_t(:, is_oce, k) = 0.
207 elsewhere
208 where (pctsrf_t(:, is_sic, k) >= 1 - zmasq)
209 pctsrf_t(:, is_sic, k) = 1. - zmasq
210 pctsrf_t(:, is_oce, k) = 0.
211 ELSEwhere
212 pctsrf_t(:, is_oce, k) = 1. - zmasq - pctsrf_t(:, is_sic, k)
213 where (pctsrf_t(:, is_oce, k) < EPSFRA)
214 pctsrf_t(:, is_oce, k) = 0.
215 pctsrf_t(:, is_sic, k) = 1 - zmasq
216 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 print *, 'Bad surface fraction: pctsrf_t(', i, &
223 ', 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 print *, 'Bad surface fraction:'
229 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 call nf95_inq_varid(ncid, "longitude", varid)
240 call nf95_gw_var(ncid, varid, dlon_ini)
241 imdep = size(dlon_ini)
242
243 call nf95_inq_varid(ncid, "latitude", varid)
244 call nf95_gw_var(ncid, varid, dlat_ini)
245 jmdep = size(dlat_ini)
246
247 call nf95_inq_dimid(ncid, "time", dimid)
248 call NF95_INQuire_DIMension(ncid, dimid, nclen=lmdep)
249 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 ALLOCATE(champ(imdep, jmdep), champtime(iim, jjm + 1, lmdep))
255 IF(extrap) THEN
256 ALLOCATE(work(imdep, jmdep))
257 ENDIF
258 ALLOCATE(dlon(imdep), dlat(jmdep))
259 call NF95_INQ_VARID(ncid, 'tosbcs', varid)
260
261 DO l = 1, lmdep
262 call NF95_GET_VAR(ncid, varid, champ, start=(/1, 1, l/))
263 CALL conf_dat2d(dlon_ini, dlat_ini, dlon, dlat, champ)
264 IF (extrap) THEN
265 CALL extrapol(champ, imdep, jmdep, 999999., .TRUE., .TRUE., 2, work)
266 ENDIF
267
268 CALL inter_barxy(dlon, dlat(:jmdep -1), champ, rlonu(:iim), rlatv, &
269 champtime(:, :, l))
270 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 yder = SPLINE(tmidmonth, champtime(i, j, :))
281 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 PRINT *, "Traitement de l'albedo..."
304 call NF95_OPEN('Albedo.nc', NF90_NOWRITE, ncid)
305
306 call nf95_inq_varid(ncid, "longitude", varid)
307 call nf95_gw_var(ncid, varid, dlon_ini)
308 imdep = size(dlon_ini)
309
310 call nf95_inq_varid(ncid, "latitude", varid)
311 call nf95_gw_var(ncid, varid, dlat_ini)
312 jmdep = size(dlat_ini)
313
314 call nf95_inq_varid(ncid, "temps", varid)
315 call nf95_gw_var(ncid, varid, timeyear)
316 lmdep = size(timeyear)
317
318 ALLOCATE(champ(imdep, jmdep), champtime(iim, jjm + 1, lmdep))
319 ALLOCATE(dlon(imdep), dlat(jmdep))
320 call NF95_INQ_VARID(ncid, 'ALBEDO', varid)
321
322 DO l = 1, lmdep
323 PRINT *, "timeyear(", l, ") =", timeyear(l)
324 call NF95_GET_VAR(ncid, varid, champ, start=(/1, 1, l/))
325 CALL conf_dat2d(dlon_ini, dlat_ini, dlon, dlat, champ)
326 CALL inter_barxy(dlon, dlat(:jmdep-1), champ, rlonu(:iim), rlatv, &
327 champtime(:, :, l))
328 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 yder = SPLINE(timeyear, champtime(i, j, :))
339 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 PRINT *, 'Ecriture du fichier limit.nc...'
357
358 call NF95_CREATE("limit.nc", NF90_CLOBBER, nid)
359
360 call NF95_PUT_ATT(nid, NF90_GLOBAL, "title", &
361 "Fichier conditions aux limites")
362 call NF95_DEF_DIM(nid, "points_physiques", klon, ndim)
363 call NF95_DEF_DIM(nid, "time", NF90_UNLIMITED, ntim)
364
365 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
368 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
372 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
376 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
380 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
384 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 "Temperature superficielle de la mer")
388
389 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 call NF95_ENDDEF(nid)
402
403 DO k = 1, 360
404 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 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