Context Navigation

create_limit.F @ 266

Last change on this file since 266 was 232, checked in by lmdzadmin, 23 years ago
Merge par rapport a la branche principale LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 32.5 KB

Line
1	c $Header$
2	PROGRAM create_limit
3	USE startvar
4	USE ioipsl
5	IMPLICIT none
6	c
7	c-------------------------------------------------------------
8	C Author : L. Fairhead
9	C Date : 27/01/94
10	C Objet : Construction des fichiers de conditions aux limites
11	C pour le nouveau
12	C modele a partir de fichiers de climatologie. Les deux
13	C grilles doivent etre regulieres
14	c
15	c Modifie par z.x.li (le23mars1994)
16	c Modifie par L. Fairhead (fairhead@lmd.jussieu.fr) septembre 1999
17	c pour lecture netcdf dans LMDZ.3.3
18	c modifie par P. Braconnot pour utiliser la version sous-surfaces
19	c-------------------------------------------------------------
20	c
21	#include "dimensions.h"
22	#include "paramet.h"
23	#include "control.h"
24	#include "logic.h"
25	#include "netcdf.inc"
26	#include "comvert.h"
27	#include "comgeom2.h"
28	#include "comconst.h"
29	#include "dimphy.h"
30	#include "indicesol.h"
31	c-----------------------------------------------------------------------
32	REAL phy_nat(klon,360)
33	real phy_nat0(klon)
34	REAL phy_alb(klon,360)
35	REAL phy_sst(klon,360)
36	REAL phy_bil(klon,360)
37	REAL phy_rug(klon,360)
38	REAL phy_ice(klon)
39	CPB
40	c REAL phy_icet(klon,360)
41	c REAL phy_oce(klon,360)
42	real pctsrf_t(klon,nbsrf,360)
43	real pctsrf(klon,nbsrf)
44	REAL verif
45	c
46	REAL masque(iip1,jjp1)
47	REAL mask(iim,jjp1)
48	CPB
49	C newlmt indique l'utilisation de la sous-maille fractionnelle
50	C tandis que l'ancien codage utilise l'indicateur du sol (0,1,2,3)
51	LOGICAL newlmt, fracterre
52	PARAMETER(newlmt=.TRUE.)
53	PARAMETER(fracterre = .TRUE.)
54	CPB
55	C Declarations pour le champ de depart
56	INTEGER imdep, jmdep,lmdep
57	INTEGER ibid, jbid, tbid
58	PARAMETER (ibid = 400, ! >360 pts
59	. jbid = 200, ! >181 pts
60	. tbid = 60) ! >52 semaines
61	REAL champ(ibid*jbid)
62	REAL dlon(ibid), dlat(jbid), timecoord(tbid)
63	c
64	INTEGER ibid_msk, jbid_msk
65	PARAMETER(ibid_msk=2200,jbid_msk=1100)
66	REAL champ_msk(ibid_msk*jbid_msk)
67	REAL dlon_msk(ibid_msk), dlat_msk(jbid_msk)
68	REAL4 zbidon(ibid_mskjbid_msk)
69	C Declarations pour le champ interpole 2D
70	REAL champint(iim,jjp1)
71
72	C Declarations pour le champ interpole 3D
73	REAL champtime(iim,jjp1,tbid)
74	REAL timeyear(tbid)
75	REAL champan(iip1,jjp1,366)
76
77	C Declarations pour l'inteprolation verticale
78	REAL ax(tbid), ay(tbid)
79	REAL by
80	REAL yder(tbid)
81
82
83	INTEGER ierr
84	INTEGER dimfirst(3)
85	INTEGER dimlast(3)
86	c
87	INTEGER nid, ndim, ntim
88	INTEGER dims(2), debut(2), epais(2)
89	INTEGER id_tim
90	INTEGER id_NAT, id_SST, id_BILS, id_RUG, id_ALB
91	CPB
92	INTEGER id_FOCE, id_FSIC, id_FTER, id_FLIC
93
94	INTEGER i, j, k, l, ji
95	c declarations pour lecture glace de mer
96	INTEGER :: iml_lic, jml_lic, llm_tmp, ttm_tmp, iret
97	INTEGER :: itaul(1), fid
98	REAL :: lev(1), date, dt
99	REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_lic, lat_lic
100	REAL, ALLOCATABLE, DIMENSION(:) :: dlon_lic, dlat_lic
101	REAL, ALLOCATABLE, DIMENSION (:,:) :: fraclic
102	REAL :: flic_tmp(iip1, jjp1)
103	c Diverses variables locales
104	REAL time
105	! pour la lecture du fichier masque ocean
106	integer :: nid_o2a
107	logical :: couple = .false.
108	INTEGER :: iml_omask, jml_omask
109	REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_omask, lat_omask
110	REAL, ALLOCATABLE, DIMENSION(:) :: dlon_omask, dlat_omask
111	REAL, ALLOCATABLE, DIMENSION (:,:) :: ocemask, ocetmp
112	real, dimension(klon) :: ocemask_fi
113
114
115	INTEGER longcles
116	PARAMETER ( longcles = 20 )
117	REAL clesphy0 ( longcles )
118	#include "serre.h"
119	INTEGER ncid,varid,ndimid(4),dimid
120	character*30 namedim
121	CHARACTER*80 :: varname
122
123	c initialisations:
124	OPEN (8,file='run.def',form='formatted')
125	CALL defrun_new(8,.TRUE.,clesphy0)
126	CLOSE(8)
127
128	pi = 4. * ATAN(1.)
129	rad = 6 371 229.
130	omeg = 4.* ASIN(1.)/(24.*3600.)
131	g = 9.8
132	daysec = 86400.
133	kappa = 0.2857143
134	cpp = 1004.70885
135	dtvr = daysec/FLOAT(day_step)
136
137	c
138	ccc CALL iniconst ( non indispensable )
139
140	CALL inigeom
141	c
142	c
143	C Traitement du relief au sol
144	c
145	write(,) 'Fabrication masque'
146
147	varname = 'masque'
148	masque(:,:) = 0.0
149	CALL startget(varname, iip1, jjp1, rlonv, rlatu, masque, 0.0)
150	pctsrf=0.
151	varname = 'zmasq'
152	zmasq(:) = 0.
153	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zmasq,0.0)
154	WHERE (zmasq(1 : klon) .LT. EPSFRA)
155	zmasq(1 : klon) = 0.
156	END WHERE
157	WHERE (1 - zmasq(1 : klon) .LT. EPSFRA)
158	zmasq(1 : klon) = 1.
159	END WHERE
160	! WRITE(,)zmasq
161
162	IF ( fracterre ) THEN
163	DO i = 1, iim
164	masque(i,1) = masque(i,1)
165	masque(i,jjp1) = masque(i,jjp1)
166	END DO
167	ELSE
168	DO i = 1, iim
169	masque(i,1) = FLOAT(NINT(masque(i,1)))
170	masque(i,jjp1) = FLOAT(NINT(masque(i,jjp1)))
171	END DO
172	ENDIF
173	c$$$ DO i = 1, iim
174	c$$$ DO j = 1, jjp1
175	c$$$ mask(i,j) = masque(i,j)
176	c$$$ ENDDO
177	c$$$ ENDDO
178	c$$$ CALL gr_dyn_fi(1, iip1, jjp1, klon, masque, phy_nat0)
179	phy_nat0(1:klon) = zmasq(1:klon)
180	mask = 0.
181	DO j = 1, jjp1
182	DO i = 1, iim
183	IF ( masque(i,j) .GE. EPSFRA) mask (i,j) = 1
184	END DO
185	END DO
186	C
187	C En cas de simulation couplee, lecture du masque ocean issu du modele ocean
188	C utilise pour calculer les poids et pour assurer l'adequation entre les
189	C fractions d'ocean vu par l'atmosphere et l'ocean
190	C
191
192	write(,)'Essai de lecture masque ocean'
193	iret = nf_open("o2a.nc", NF_NOWRITE, nid_o2a)
194	if (iret .ne. 0) then
195	write(,)'ATTENTION!! pas de fichier o2a.nc trouve'
196	write(,)'Run force'
197	else
198	couple = .true.
199	iret = nf_close(nid_o2a)
200	call flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp
201	$ , nid_o2a)
202	if (iml_omask /= iim .or. jml_omask /= jjp1) then
203	write(,)'Dimensions non compatibles pour masque ocean'
204	write(,)'iim = ',iim,' iml_omask = ',iml_omask
205	write(,)'jjp1 = ',jjp1,' jml_omask = ',jml_omask
206	stop
207	endif
208	ALLOCATE(lat_omask(iml_omask, jml_omask), stat=iret)
209	ALLOCATE(lon_omask(iml_omask, jml_omask), stat=iret)
210	ALLOCATE(dlon_omask(iml_omask), stat=iret)
211	ALLOCATE(dlat_omask(jml_omask), stat=iret)
212	ALLOCATE(ocemask(iml_omask, jml_omask), stat=iret)
213	ALLOCATE(ocetmp(iml_omask, jml_omask), stat=iret)
214	CALL flinopen("o2a.nc", .FALSE., iml_omask, jml_omask, llm_tmp
215	$ , lon_omask, lat_omask, lev, ttm_tmp, itaul, date, dt, fid)
216	CALL flinget(fid, 'OceMask', iml_omask, jml_omask, llm_tmp,
217	$ ttm_tmp, 1, 1, ocetmp)
218	CALL flinclo(fid)
219	dlon_omask(1 : iml_omask) = lon_omask(1 : iml_omask, 1)
220	dlat_omask(1 : jml_omask) = lat_omask(1 , 1 : jml_omask)
221	ocemask = ocetmp
222	if (dlat_omask(1) < dlat_omask(jml_omask)) then
223	do j = 1, jml_omask
224	ocemask(:,j) = ocetmp(:,jml_omask-j+1)
225	enddo
226	endif
227	C
228	C passage masque ocean a la grille physique
229	C
230	ocemask_fi(1) = ocemask(1,1)
231	do j = 2, jjm
232	do i = 1, iim
233	ocemask_fi((j-2)*iim + i + 1) = ocemask(i,j)
234	enddo
235	enddo
236	ocemask_fi(klon) = ocemask(1,jjp1)
237	zmasq = 1. - ocemask_fi
238	endif
239
240
241	C
242	C lecture du fichier glace de terre pour fixer la fraction de terre
243	C et de glace de terre
244	C
245	CALL flininfo("landiceref.nc", iml_lic, jml_lic,llm_tmp, ttm_tmp
246	$ , fid)
247	ALLOCATE(lat_lic(iml_lic, jml_lic), stat=iret)
248	ALLOCATE(lon_lic(iml_lic, jml_lic), stat=iret)
249	ALLOCATE(dlon_lic(iml_lic), stat=iret)
250	ALLOCATE(dlat_lic(jml_lic), stat=iret)
251	ALLOCATE(fraclic(iml_lic, jml_lic), stat=iret)
252	CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp
253	$ , lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid)
254	CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp
255	$ , 1, 1, fraclic)
256	CALL flinclo(fid)
257	C
258	C interpolation sur la grille T du modele
259	C
260	WRITE(,) 'dimensions de landice iml_lic, jml_lic : ',
261	$ iml_lic, jml_lic
262	c
263	C sil les coordonnees sont en degres, on les transforme
264	C
265	IF( MAXVAL( lon_lic(:,:) ) .GT. 2.0 * asin(1.0) ) THEN
266	lon_lic(:,:) = lon_lic(:,:) * 2.0* ASIN(1.0) / 180.
267	ENDIF
268	IF( maxval( lat_lic(:,:) ) .GT. 2.0 * asin(1.0)) THEN
269	lat_lic(:,:) = lat_lic(:,:) * 2.0 * asin(1.0) / 180.
270	ENDIF
271
272	dlon_lic(1 : iml_lic) = lon_lic(1 : iml_lic, 1)
273	dlat_lic(1 : jml_lic) = lat_lic(1 , 1 : jml_lic)
274	C
275	CALL grille_m(iml_lic, jml_lic, dlon_lic, dlat_lic, fraclic
276	$ ,iim, jjp1,
277	$ rlonv, rlatu, flic_tmp(1 : iim, 1 : jjp1))
278	c$$$ flic_tmp(1 : iim, 1 : jjp1) = champint(1: iim, 1 : jjp1)
279	flic_tmp(iip1, 1 : jjp1) = flic_tmp(1 , 1 : jjp1)
280	C
281	C passage sur la grille physique
282	C
283	CALL gr_dyn_fi(1, iip1, jjp1, klon, flic_tmp,
284	$ pctsrf(1:klon, is_lic))
285	C adequation avec le maque terre/mer
286	WHERE (pctsrf(1 : klon, is_lic) .LT. EPSFRA )
287	pctsrf(1 : klon, is_lic) = 0.
288	END WHERE
289	WHERE (zmasq( 1 : klon) .LT. EPSFRA)
290	pctsrf(1 : klon, is_lic) = 0.
291	END WHERE
292	pctsrf(1 : klon, is_ter) = zmasq(1 : klon)
293	DO ji = 1, klon
294	IF (zmasq(ji) .GT. EPSFRA) THEN
295	IF ( pctsrf(ji, is_lic) .GE. zmasq(ji)) THEN
296	pctsrf(ji, is_lic) = zmasq(ji)
297	pctsrf(ji, is_ter) = 0.
298	ELSE
299	pctsrf(ji,is_ter) = zmasq(ji) - pctsrf(ji, is_lic)
300	IF (pctsrf(ji,is_ter) .LT. EPSFRA) THEN
301	pctsrf(ji,is_ter) = 0.
302	pctsrf(ji, is_lic) = zmasq(ji)
303	ENDIF
304	ENDIF
305	ENDIF
306	END DO
307	c
308	c
309	C Traitement de la rugosite
310	c
311	PRINT*, 'Traitement de la rugosite'
312	ierr = NF_OPEN('Rugos.nc', NF_NOWRITE, ncid)
313	if (ierr.ne.0) then
314	print *, NF_STRERROR(ierr)
315	STOP
316	ENDIF
317
318	ierr = NF_INQ_VARID(ncid,'RUGOS',varid)
319	if (ierr.ne.0) then
320	print *, NF_STRERROR(ierr)
321	STOP
322	ENDIF
323	ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
324	if (ierr.ne.0) then
325	print *, NF_STRERROR(ierr)
326	STOP
327	ENDIF
328	ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
329	if (ierr.ne.0) then
330	print *, NF_STRERROR(ierr)
331	STOP
332	ENDIF
333	print*,'variable ', namedim, 'dimension ', imdep
334	ierr = NF_INQ_VARID(ncid,namedim,dimid)
335	if (ierr.ne.0) then
336	print *, NF_STRERROR(ierr)
337	STOP
338	ENDIF
339	#ifdef NC_DOUBLE
340	ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon)
341	#else
342	ierr = NF_GET_VAR_REAL(ncid,dimid,dlon)
343	#endif
344	if (ierr.ne.0) then
345	print *, NF_STRERROR(ierr)
346	STOP
347	ENDIF
348	ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
349	if (ierr.ne.0) then
350	print *, NF_STRERROR(ierr)
351	STOP
352	ENDIF
353	print*,'variable ', namedim, 'dimension ', jmdep
354	ierr = NF_INQ_VARID(ncid,namedim,dimid)
355	if (ierr.ne.0) then
356	print *, NF_STRERROR(ierr)
357	STOP
358	ENDIF
359	#ifdef NC_DOUBLE
360	ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat)
361	#else
362	ierr = NF_GET_VAR_REAL(ncid,dimid,dlat)
363	#endif
364	if (ierr.ne.0) then
365	print *, NF_STRERROR(ierr)
366	STOP
367	ENDIF
368	ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep)
369	if (ierr.ne.0) then
370	print *, NF_STRERROR(ierr)
371	STOP
372	ENDIF
373	print*,'variable ', namedim, 'dimension ', lmdep
374	ierr = NF_INQ_VARID(ncid,namedim,dimid)
375	if (ierr.ne.0) then
376	print *, NF_STRERROR(ierr)
377	STOP
378	ENDIF
379	#ifdef NC_DOUBLE
380	ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
381	#else
382	ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
383	#endif
384	if (ierr.ne.0) then
385	print *, NF_STRERROR(ierr)
386	STOP
387	ENDIF
388	c
389	DO l = 1, lmdep
390	dimfirst(1) = 1
391	dimfirst(2) = 1
392	dimfirst(3) = l
393	c
394	dimlast(1) = imdep
395	dimlast(2) = jmdep
396	dimlast(3) = 1
397	c
398	PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
399	print*,dimfirst,dimlast
400	#ifdef NC_DOUBLE
401	ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
402	#else
403	ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
404	#endif
405	if (ierr.ne.0) then
406	print *, NF_STRERROR(ierr)
407	STOP
408	ENDIF
409
410	CALL rugosite(imdep, jmdep, dlon, dlat, champ,
411	. iim, jjp1, rlonv, rlatu, champint, mask)
412	DO j = 1,jjp1
413	DO i = 1, iim
414	champtime (i,j,l) = champint(i,j)
415	ENDDO
416	ENDDO
417	ENDDO
418	c write(70,*)champtime
419	c
420	DO l = 1, lmdep
421	timeyear(l) = timecoord(l)
422	ENDDO
423
424	PRINT 222, timeyear
425	222 FORMAT(2x,' Time year ',10f6.1)
426	c
427
428	PRINT*, 'Interpolation temporelle dans l annee'
429
430
431	DO j = 1, jjp1
432	DO i = 1, iim
433	DO l = 1, lmdep
434	ax(l) = timeyear(l)
435	ay(l) = champtime (i,j,l)
436	ENDDO
437	CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder)
438	DO k = 1, 360
439	time = FLOAT(k-1)
440	CALL SPLINT(ax,ay,yder,lmdep,time,by)
441	champan(i,j,k) = by
442	ENDDO
443	ENDDO
444	ENDDO
445	DO k = 1, 360
446	DO j = 1, jjp1
447	champan(iip1,j,k) = champan(1,j,k)
448	ENDDO
449	ENDDO
450	c
451	DO k = 1, 360
452	CALL gr_dyn_fi(1,iip1,jjp1,klon,champan(1,1,k), phy_rug(1,k))
453	ENDDO
454	c
455	ierr = NF_CLOSE(ncid)
456	c
457	c
458	C Traitement de la glace oceanique
459	c
460	PRINT*, 'Traitement de la glace oceanique'
461	ierr = NF_OPEN('AMIP.nc', NF_NOWRITE, ncid)
462	if (ierr.ne.0) then
463	print *, NF_STRERROR(ierr)
464	STOP
465	ENDIF
466
467	ierr = NF_INQ_VARID(ncid,'SEA_ICE',varid)
468	if (ierr.ne.0) then
469	print *, NF_STRERROR(ierr)
470	STOP
471	ENDIF
472	ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
473	if (ierr.ne.0) then
474	print *, NF_STRERROR(ierr)
475	STOP
476	ENDIF
477	ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
478	if (ierr.ne.0) then
479	print *, NF_STRERROR(ierr)
480	STOP
481	ENDIF
482	print*,'variable ', namedim, 'dimension ', imdep
483	ierr = NF_INQ_VARID(ncid,namedim,dimid)
484	if (ierr.ne.0) then
485	print *, NF_STRERROR(ierr)
486	STOP
487	ENDIF
488	#ifdef NC_DOUBLE
489	ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon)
490	#else
491	ierr = NF_GET_VAR_REAL(ncid,dimid,dlon)
492	#endif
493	if (ierr.ne.0) then
494	print *, NF_STRERROR(ierr)
495	STOP
496	ENDIF
497	ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
498	if (ierr.ne.0) then
499	print *, NF_STRERROR(ierr)
500	STOP
501	ENDIF
502	print*,'variable ', namedim, jmdep
503	ierr = NF_INQ_VARID(ncid,namedim,dimid)
504	if (ierr.ne.0) then
505	print *, NF_STRERROR(ierr)
506	STOP
507	ENDIF
508	#ifdef NC_DOUBLE
509	ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat)
510	#else
511	ierr = NF_GET_VAR_REAL(ncid,dimid,dlat)
512	#endif
513	if (ierr.ne.0) then
514	print *, NF_STRERROR(ierr)
515	STOP
516	ENDIF
517	ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep)
518	if (ierr.ne.0) then
519	print *, NF_STRERROR(ierr)
520	STOP
521	ENDIF
522	print*,'variable ', namedim, lmdep
523	ierr = NF_INQ_VARID(ncid,namedim,dimid)
524	if (ierr.ne.0) then
525	print *, NF_STRERROR(ierr)
526	STOP
527	ENDIF
528	#ifdef NC_DOUBLE
529	ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
530	#else
531	ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
532	#endif
533	if (ierr.ne.0) then
534	print *, NF_STRERROR(ierr)
535	STOP
536	ENDIF
537	c
538	DO l = 1, lmdep
539	dimfirst(1) = 1
540	dimfirst(2) = 1
541	dimfirst(3) = l
542	c
543	dimlast(1) = imdep
544	dimlast(2) = jmdep
545	dimlast(3) = 1
546	c
547	PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
548	#ifdef NC_DOUBLE
549	ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
550	#else
551	ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
552	#endif
553	if (ierr.ne.0) then
554	print *, NF_STRERROR(ierr)
555	STOP
556	ENDIF
557
558	CALL sea_ice(imdep, jmdep, dlon, dlat, champ,
559	. iim, jjp1, rlonv, rlatu, champint)
560	DO j = 1,jjp1
561	DO i = 1, iim
562	champtime (i,j,l) = champint(i,j)
563	ENDDO
564	ENDDO
565	ENDDO
566	c
567	DO l = 1, lmdep
568	timeyear(l) = timecoord(l)
569	ENDDO
570	PRINT 222, timeyear
571	c
572	PRINT*, 'Interpolation temporelle'
573	DO j = 1, jjp1
574	DO i = 1, iim
575	DO l = 1, lmdep
576	ax(l) = timeyear(l)
577	ay(l) = champtime (i,j,l)
578	ENDDO
579	CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder)
580	DO k = 1, 360
581	time = FLOAT(k-1)
582	CALL SPLINT(ax,ay,yder,lmdep,time,by)
583	champan(i,j,k) = by
584	ENDDO
585	ENDDO
586	ENDDO
587	DO k = 1, 360
588	DO j = 1, jjp1
589	champan(iip1, j, k) = champan(1, j, k)
590	ENDDO
591	ENDDO
592	c
593	c WRITE(,) 'phy_nat'
594	c WRITE(*,'(72f4.1)') phy_nat0(1:klon)
595	c
596	DO k = 1, 360
597	CALL gr_dyn_fi(1, iip1, jjp1, klon,
598	. champan(1,1,k), phy_ice)
599	IF ( newlmt) THEN
600
601	CPB en attendant de mettre fraction de terre
602	c
603	WHERE(phy_ice(1:klon) .GE. 1.) phy_ice(1 : klon) = 1.
604	WHERE(phy_ice(1:klon) .LT. EPSFRA) phy_ice(1 : klon) = 0.
605	c
606	IF (fracterre ) THEN
607	c WRITE(,) 'passe dans cas fracterre'
608	pctsrf_t(:,is_ter,k) = pctsrf(:,is_ter)
609	pctsrf_t(:,is_lic,k) = pctsrf(:,is_lic)
610	pctsrf_t(1:klon,is_sic,k) = phy_ice(1:klon)
611	$ - pctsrf_t(1:klon,is_lic,k)
612	c§§ Il y a des cas ou il y a de la glace dans landiceref et pas dans AMIP
613	WHERE (pctsrf_t(1:klon,is_sic,k) .LE. 0)
614	pctsrf_t(1:klon,is_sic,k) = 0.
615	END WHERE
616	WHERE( 1. - zmasq(1:klon) .LT. EPSFRA)
617	pctsrf_t(1:klon,is_sic,k) = 0.
618	pctsrf_t(1:klon,is_oce,k) = 0.
619	END WHERE
620	DO i = 1, klon
621	c$$ pctsrf_t(i,is_sic,k) = (1. - pctsrf_t(i,is_lic,k) -
622	c$$ . pctsrf_t(i,is_ter,k)) * phy_ice(i)
623	c$$ pctsrf_t(i,is_oce,k) = 1. - pctsrf_t(i,is_lic,k) -
624	c$$ . pctsrf_t(i,is_ter,k) - pctsrf_t(i,is_sic,k)
625	IF ( 1. - zmasq(i) .GT. EPSFRA) THEN
626	IF ( pctsrf_t(i,is_sic,k) .GE. 1 - zmasq(i)) THEN
627	pctsrf_t(i,is_sic,k) = 1 - zmasq(i)
628	pctsrf_t(i,is_oce,k) = 0.
629	ELSE
630	pctsrf_t(i,is_oce,k) = 1 - zmasq(i)
631	$ - pctsrf_t(i,is_sic,k)
632	IF (pctsrf_t(i,is_oce,k) .LT. EPSFRA) THEN
633	pctsrf_t(i,is_oce,k) = 0.
634	pctsrf_t(i,is_sic,k) = 1 - zmasq(i)
635	ENDIF
636	ENDIF
637	ENDIF
638	if (pctsrf_t(i,is_oce,k) .lt. 0.) then
639	WRITE(,) 'pb sous maille au point : i,k '
640	$ , i,k,pctsrf_t(:,is_oce,k)
641	ENDIF
642	IF ( abs( pctsrf_t(i, is_ter,k) + pctsrf_t(i, is_lic,k) +
643	$ pctsrf_t(i, is_oce,k) + pctsrf_t(i, is_sic,k) - 1.)
644	$ .GT. EPSFRA) THEN
645	WRITE(,) 'physiq : pb sous surface au point ', i,
646	$ pctsrf_t(i, 1 : nbsrf,k), phy_ice(i)
647	ENDIF
648	END DO
649	ELSE
650	DO i = 1, klon
651	pctsrf_t(i,is_ter,k) = pctsrf(i,is_ter)
652	IF (NINT(pctsrf(i,is_ter)).EQ.1 ) THEN
653	pctsrf_t(i,is_sic,k) = 0.
654	pctsrf_t(i,is_oce,k) = 0.
655	IF(phy_ice(i) .GE. 1.e-5) THEN
656	pctsrf_t(i,is_lic,k) = phy_ice(i)
657	pctsrf_t(i,is_ter,k) = pctsrf_t(i,is_ter,k)
658	. - pctsrf_t(i,is_lic,k)
659	ELSE
660	pctsrf_t(i,is_lic,k) = 0.
661	ENDIF
662	ELSE
663	pctsrf_t(i,is_lic,k) = 0.
664	IF(phy_ice(i) .GE. 1.e-5) THEN
665	pctsrf_t(i,is_ter,k) = 0.
666	pctsrf_t(i,is_sic,k) = phy_ice(i)
667	pctsrf_t(i,is_oce,k) = 1. - pctsrf_t(i,is_sic,k)
668	ELSE
669	pctsrf_t(i,is_sic,k) = 0.
670	pctsrf_t(i,is_oce,k) = 1.
671	ENDIF
672	ENDIF
673	verif = pctsrf_t(i,is_ter,k) +
674	. pctsrf_t(i,is_oce,k) +
675	. pctsrf_t(i,is_sic,k) +
676	. pctsrf_t(i,is_lic,k)
677	IF ( verif .LT. 1. - 1.e-5 .OR.
678	$ verif .GT. 1 + 1.e-5) THEN
679	WRITE(,) 'pb sous maille au point : i,k,verif '
680	$ , i,k,verif
681	ENDIF
682	END DO
683	ENDIF
684	ELSE
685	DO i = 1, klon
686	phy_nat(i,k) = phy_nat0(i)
687	IF ( (phy_ice(i) - 0.5).GE.1.e-5 ) THEN
688	IF (NINT(phy_nat0(i)).EQ.0) THEN
689	phy_nat(i,k) = 3.0
690	ELSE
691	phy_nat(i,k) = 2.0
692	ENDIF
693	ENDIF
694	END DO
695	ENDIF
696	ENDDO
697	c
698	ierr = NF_CLOSE(ncid)
699	c
700	c
701	C Traitement de la sst
702	c
703	PRINT*, 'Traitement de la sst'
704	ierr = NF_OPEN('AMIP.nc', NF_NOWRITE, ncid)
705	if (ierr.ne.0) then
706	print *, NF_STRERROR(ierr)
707	STOP
708	ENDIF
709
710	ierr = NF_INQ_VARID(ncid,'SST',varid)
711	if (ierr.ne.0) then
712	print *, NF_STRERROR(ierr)
713	STOP
714	ENDIF
715	ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
716	if (ierr.ne.0) then
717	print *, NF_STRERROR(ierr)
718	STOP
719	ENDIF
720	ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
721	if (ierr.ne.0) then
722	print *, NF_STRERROR(ierr)
723	STOP
724	ENDIF
725	print*,'variable ', namedim,'dimension ', imdep
726	ierr = NF_INQ_VARID(ncid,namedim,dimid)
727	if (ierr.ne.0) then
728	print *, NF_STRERROR(ierr)
729	STOP
730	ENDIF
731	#ifdef NC_DOUBLE
732	ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon)
733	#else
734	ierr = NF_GET_VAR_REAL(ncid,dimid,dlon)
735	#endif
736
737	if (ierr.ne.0) then
738	print *, NF_STRERROR(ierr)
739	STOP
740	ENDIF
741	ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
742	if (ierr.ne.0) then
743	print *, NF_STRERROR(ierr)
744	STOP
745	ENDIF
746	print*,'variable ', namedim, 'dimension ', jmdep
747	ierr = NF_INQ_VARID(ncid,namedim,dimid)
748	if (ierr.ne.0) then
749	print *, NF_STRERROR(ierr)
750	STOP
751	ENDIF
752	#ifdef NC_DOUBLE
753	ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat)
754	#else
755	ierr = NF_GET_VAR_REAL(ncid,dimid,dlat)
756	#endif
757	if (ierr.ne.0) then
758	print *, NF_STRERROR(ierr)
759	STOP
760	ENDIF
761	ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep)
762	if (ierr.ne.0) then
763	print *, NF_STRERROR(ierr)
764	STOP
765	ENDIF
766	print*,'variable ', namedim, 'dimension ', lmdep
767	ierr = NF_INQ_VARID(ncid,namedim,dimid)
768	if (ierr.ne.0) then
769	print *, NF_STRERROR(ierr)
770	STOP
771	ENDIF
772	#ifdef NC_DOUBLE
773	ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
774	#else
775	ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
776	#endif
777	if (ierr.ne.0) then
778	print *, NF_STRERROR(ierr)
779	STOP
780	ENDIF
781	c
782	DO l = 1, lmdep
783	dimfirst(1) = 1
784	dimfirst(2) = 1
785	dimfirst(3) = l
786	c
787	dimlast(1) = imdep
788	dimlast(2) = jmdep
789	dimlast(3) = 1
790	c
791	PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
792	#ifdef NC_DOUBLE
793	ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
794	#else
795	ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
796	#endif
797	if (ierr.ne.0) then
798	print *, NF_STRERROR(ierr)
799	STOP
800	ENDIF
801	CALL grille_m(imdep, jmdep, dlon, dlat, champ,
802	. iim, jjp1, rlonv, rlatu, champint)
803
804	DO j = 1,jjp1
805	DO i = 1, iim
806	champtime (i,j,l) = champint(i,j)
807	ENDDO
808	ENDDO
809	ENDDO
810	c
811	DO l = 1, lmdep
812	timeyear(l) = timecoord(l)
813	ENDDO
814	print 222, timeyear
815	c
816	C interpolation temporelle
817	DO j = 1, jjp1
818	DO i = 1, iim
819	DO l = 1, lmdep
820	ax(l) = timeyear(l)
821	ay(l) = champtime (i,j,l)
822	ENDDO
823	CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder)
824	DO k = 1, 360
825	time = FLOAT(k-1)
826	CALL SPLINT(ax,ay,yder,lmdep,time,by)
827	champan(i,j,k) = by
828	ENDDO
829	ENDDO
830	ENDDO
831	DO k = 1, 360
832	DO j = 1, jjp1
833	champan(iip1,j,k) = champan(1,j,k)
834	ENDDO
835	ENDDO
836	c
837	DO k = 1, 360
838	CALL gr_dyn_fi(1, iip1, jjp1, klon,
839	. champan(1,1,k), phy_sst(1,k))
840	ENDDO
841	c
842	WHERE(phy_sst .LT. 271.35) phy_sst = 271.35
843	ierr = NF_CLOSE(ncid)
844	c
845	c
846	C Traitement de l'albedo
847	c
848	PRINT*, 'Traitement de l albedo'
849	ierr = NF_OPEN('Albedo.nc', NF_NOWRITE, ncid)
850	if (ierr.ne.0) then
851	print *, NF_STRERROR(ierr)
852	STOP
853	ENDIF
854	ierr = NF_INQ_VARID(ncid,'ALBEDO',varid)
855	if (ierr.ne.0) then
856	print *, NF_STRERROR(ierr)
857	STOP
858	ENDIF
859	ierr = NF_INQ_VARDIMID (ncid,varid,ndimid)
860	if (ierr.ne.0) then
861	print *, NF_STRERROR(ierr)
862	STOP
863	ENDIF
864	ierr = NF_INQ_DIM(ncid,ndimid(1), namedim, imdep)
865	if (ierr.ne.0) then
866	print *, NF_STRERROR(ierr)
867	STOP
868	ENDIF
869	print*,'variable ', namedim, 'dimension ', imdep
870	ierr = NF_INQ_VARID(ncid,namedim,dimid)
871	if (ierr.ne.0) then
872	print *, NF_STRERROR(ierr)
873	STOP
874	ENDIF
875	#ifdef NC_DOUBLE
876	ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlon)
877	#else
878	ierr = NF_GET_VAR_REAL(ncid,dimid,dlon)
879	#endif
880	if (ierr.ne.0) then
881	print *, NF_STRERROR(ierr)
882	STOP
883	ENDIF
884	ierr = NF_INQ_DIM(ncid,ndimid(2), namedim, jmdep)
885	if (ierr.ne.0) then
886	print *, NF_STRERROR(ierr)
887	STOP
888	ENDIF
889	print*,'variable ', namedim, 'dimension ', jmdep
890	ierr = NF_INQ_VARID(ncid,namedim,dimid)
891	if (ierr.ne.0) then
892	print *, NF_STRERROR(ierr)
893	STOP
894	ENDIF
895	#ifdef NC_DOUBLE
896	ierr = NF_GET_VAR_DOUBLE(ncid,dimid,dlat)
897	#else
898	ierr = NF_GET_VAR_REAL(ncid,dimid,dlat)
899	#endif
900	if (ierr.ne.0) then
901	print *, NF_STRERROR(ierr)
902	STOP
903	ENDIF
904	ierr = NF_INQ_DIM(ncid,ndimid(3), namedim, lmdep)
905	if (ierr.ne.0) then
906	print *, NF_STRERROR(ierr)
907	STOP
908	ENDIF
909	print*,'variable ', namedim, 'dimension ', lmdep
910	ierr = NF_INQ_VARID(ncid,namedim,dimid)
911	if (ierr.ne.0) then
912	print *, NF_STRERROR(ierr)
913	STOP
914	ENDIF
915	#ifdef NC_DOUBLE
916	ierr = NF_GET_VAR_DOUBLE(ncid,dimid,timecoord)
917	#else
918	ierr = NF_GET_VAR_REAL(ncid,dimid,timecoord)
919	#endif
920	if (ierr.ne.0) then
921	print *, NF_STRERROR(ierr)
922	STOP
923	ENDIF
924	c
925	DO l = 1, lmdep
926	dimfirst(1) = 1
927	dimfirst(2) = 1
928	dimfirst(3) = l
929	c
930	dimlast(1) = imdep
931	dimlast(2) = jmdep
932	dimlast(3) = 1
933	c
934	PRINT*,'Lecture temporelle et int. horizontale ',l,timecoord(l)
935	#ifdef NC_DOUBLE
936	ierr = NF_GET_VARA_DOUBLE(ncid,varid,dimfirst,dimlast,champ)
937	#else
938	ierr = NF_GET_VARA_REAL(ncid,varid,dimfirst,dimlast,champ)
939	#endif
940	if (ierr.ne.0) then
941	print *, NF_STRERROR(ierr)
942	STOP
943	ENDIF
944	CALL grille_m(imdep, jmdep, dlon, dlat, champ,
945	. iim, jjp1, rlonv, rlatu, champint)
946	c
947	DO j = 1,jjp1
948	DO i = 1, iim
949	champtime (i, j, l) = champint(i, j)
950	ENDDO
951	ENDDO
952	ENDDO
953	c
954	DO l = 1, lmdep
955	timeyear(l) = timecoord(l)
956	ENDDO
957	print 222, timeyear
958	c
959	C interpolation temporelle
960	DO j = 1, jjp1
961	DO i = 1, iim
962	DO l = 1, lmdep
963	ax(l) = timeyear(l)
964	ay(l) = champtime (i, j, l)
965	ENDDO
966	CALL SPLINE(ax,ay,lmdep,1.e30,1.e30,yder)
967	DO k = 1, 360
968	time = FLOAT(k-1)
969	CALL SPLINT(ax,ay,yder,lmdep,time,by)
970	champan(i,j,k) = by
971	ENDDO
972	ENDDO
973	ENDDO
974	DO k = 1, 360
975	DO j = 1, jjp1
976	champan(iip1, j, k) = champan(1, j, k)
977	ENDDO
978	ENDDO
979	c
980	DO k = 1, 360
981	CALL gr_dyn_fi(1, iip1, jjp1, klon,
982	. champan(1,1,k), phy_alb(1,k))
983	ENDDO
984	c
985	ierr = NF_CLOSE(ncid)
986	c
987	c
988	DO k = 1, 360
989	DO i = 1, klon
990	phy_bil(i,k) = 0.0
991	ENDDO
992	ENDDO
993	c
994	PRINT*, 'Ecriture du fichier limit'
995	c
996	ierr = NF_CREATE ("limit.nc", NF_CLOBBER, nid)
997	c
998	ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 30,
999	. "Fichier conditions aux limites")
1000	ierr = NF_DEF_DIM (nid, "points_physiques", klon, ndim)
1001	ierr = NF_DEF_DIM (nid, "time", NF_UNLIMITED, ntim)
1002	c
1003	dims(1) = ndim
1004	dims(2) = ntim
1005	c
1006	ierr = NF_DEF_VAR (nid, "TEMPS", NF_FLOAT, 1,ntim, id_tim)
1007	ierr = NF_PUT_ATT_TEXT (nid, id_tim, "title", 17,
1008	. "Jour dans l annee")
1009	IF (newlmt) THEN
1010	c
1011	ierr = NF_DEF_VAR (nid, "FOCE", NF_FLOAT, 2,dims, id_FOCE)
1012	ierr = NF_PUT_ATT_TEXT (nid, id_FOCE, "title", 14,
1013	. "Fraction ocean")
1014	c
1015	ierr = NF_DEF_VAR (nid, "FSIC", NF_FLOAT, 2,dims, id_FSIC)
1016	ierr = NF_PUT_ATT_TEXT (nid, id_FSIC, "title", 21,
1017	. "Fraction glace de mer")
1018	c
1019	ierr = NF_DEF_VAR (nid, "FTER", NF_FLOAT, 2,dims, id_FTER)
1020	ierr = NF_PUT_ATT_TEXT (nid, id_FTER, "title", 14,
1021	. "Fraction terre")
1022	c
1023	ierr = NF_DEF_VAR (nid, "FLIC", NF_FLOAT, 2,dims, id_FLIC)
1024	ierr = NF_PUT_ATT_TEXT (nid, id_FLIC, "title", 17,
1025	. "Fraction land ice")
1026	c
1027	ELSE
1028	ierr = NF_DEF_VAR (nid, "NAT", NF_FLOAT, 2,dims, id_NAT)
1029	ierr = NF_PUT_ATT_TEXT (nid, id_NAT, "title", 23,
1030	. "Nature du sol (0,1,2,3)")
1031	ENDIF
1032	C
1033	ierr = NF_DEF_VAR (nid, "SST", NF_FLOAT, 2,dims, id_SST)
1034	ierr = NF_PUT_ATT_TEXT (nid, id_SST, "title", 35,
1035	. "Temperature superficielle de la mer")
1036	ierr = NF_DEF_VAR (nid, "BILS", NF_FLOAT, 2,dims, id_BILS)
1037	ierr = NF_PUT_ATT_TEXT (nid, id_BILS, "title", 32,
1038	. "Reference flux de chaleur au sol")
1039	ierr = NF_DEF_VAR (nid, "ALB", NF_FLOAT, 2,dims, id_ALB)
1040	ierr = NF_PUT_ATT_TEXT (nid, id_ALB, "title", 19,
1041	. "Albedo a la surface")
1042	ierr = NF_DEF_VAR (nid, "RUG", NF_FLOAT, 2,dims, id_RUG)
1043	ierr = NF_PUT_ATT_TEXT (nid, id_RUG, "title", 8,
1044	. "Rugosite")
1045	c
1046	ierr = NF_ENDDEF(nid)
1047	c
1048	DO k = 1, 360
1049	c
1050	debut(1) = 1
1051	debut(2) = k
1052	epais(1) = klon
1053	epais(2) = 1
1054	c
1055	#ifdef NC_DOUBLE
1056	ierr = NF_PUT_VAR1_DOUBLE (nid,id_tim,k,DBLE(k))
1057	c
1058	IF (newlmt ) THEN
1059	ierr = NF_PUT_VARA_DOUBLE (nid,id_FOCE,debut,epais
1060	$ ,pctsrf_t(1,is_oce,k))
1061	ierr = NF_PUT_VARA_DOUBLE (nid,id_FSIC,debut,epais
1062	$ ,pctsrf_t(1,is_sic,k))
1063	ierr = NF_PUT_VARA_DOUBLE (nid,id_FTER,debut,epais
1064	$ ,pctsrf_t(1,is_ter,k))
1065	ierr = NF_PUT_VARA_DOUBLE (nid,id_FLIC,debut,epais
1066	$ ,pctsrf_t(1,is_lic,k))
1067	ELSE
1068	ierr = NF_PUT_VARA_DOUBLE (nid,id_NAT,debut,epais
1069	$ ,phy_nat(1,k))
1070	ENDIF
1071	c
1072	ierr = NF_PUT_VARA_DOUBLE (nid,id_SST,debut,epais,phy_sst(1,k))
1073	ierr = NF_PUT_VARA_DOUBLE (nid,id_BILS,debut,epais,phy_bil(1,k))
1074	ierr = NF_PUT_VARA_DOUBLE (nid,id_ALB,debut,epais,phy_alb(1,k))
1075	ierr = NF_PUT_VARA_DOUBLE (nid,id_RUG,debut,epais,phy_rug(1,k))
1076	#else
1077	ierr = NF_PUT_VAR1_REAL (nid,id_tim,k,FLOAT(k))
1078	IF (newlmt ) THEN
1079	ierr = NF_PUT_VARA_REAL (nid,id_FOCE,debut,epais
1080	$ ,pctsrf_t(1,is_oce,k))
1081	ierr = NF_PUT_VARA_REAL (nid,id_FSIC,debut,epais
1082	$ ,pctsrf_t(1,is_sic,k))
1083	ierr = NF_PUT_VARA_REAL (nid,id_FTER,debut,epais
1084	$ ,pctsrf_t(1,is_ter,k))
1085	ierr = NF_PUT_VARA_REAL (nid,id_FLIC,debut,epais
1086	$ ,pctsrf_t(1,is_lic,k))
1087	ELSE
1088	ierr = NF_PUT_VARA_REAL (nid,id_NAT,debut,epais
1089	$ ,phy_nat(1,k))
1090	ENDIF
1091	ierr = NF_PUT_VARA_REAL (nid,id_SST,debut,epais,phy_sst(1,k))
1092	ierr = NF_PUT_VARA_REAL (nid,id_BILS,debut,epais,phy_bil(1,k))
1093	ierr = NF_PUT_VARA_REAL (nid,id_ALB,debut,epais,phy_alb(1,k))
1094	ierr = NF_PUT_VARA_REAL (nid,id_RUG,debut,epais,phy_rug(1,k))
1095	#endif
1096	c
1097	ENDDO
1098	c
1099	ierr = NF_CLOSE(nid)
1100	c
1101	STOP
1102	END
1103

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format