Context Navigation

etat0_netcdf.F @ 724

Last change on this file since 724 was 724, checked in by Laurent Fairhead, 18 years ago
Parenthese mal placee LG LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 22.1 KB

Line
1	!
2	! $Header$
3	!
4	c
5	c
6	SUBROUTINE etat0_netcdf (interbar, masque, pctsrf)
7
8	USE startvar
9	USE ioipsl
10	!
11	IMPLICIT NONE
12	!
13	#include "netcdf.inc"
14	#include "dimensions.h"
15	#include "paramet.h"
16	!
17	!
18	! INTEGER, PARAMETER :: KIDIA=1, KFDIA=iim*(jjm-1)+2,
19	! .KLON=KFDIA-KIDIA+1,KLEV=llm
20	!
21	#include "comgeom2.h"
22	#include "comvert.h"
23	#include "comconst.h"
24	#include "indicesol.h"
25	#include "dimphy.h"
26	#include "dimsoil.h"
27	#include "temps.h"
28	!
29	LOGICAL interbar
30	REAL :: latfi(klon), lonfi(klon)
31	REAL :: orog(iip1,jjp1), rugo(iip1,jjp1), masque(iip1,jjp1),
32	. psol(iip1, jjp1), phis(iip1, jjp1)
33	REAL :: p3d(iip1, jjp1, llm+1)
34	REAL :: uvent(iip1, jjp1, llm)
35	REAL :: vvent(iip1, jjm, llm)
36	REAL :: t3d(iip1, jjp1, llm), tpot(iip1, jjp1, llm)
37	REAL :: q3d(iip1, jjp1, llm,nqmx), qsat(iip1, jjp1, llm)
38	REAL :: tsol(klon), qsol(klon), sn(klon)
39	REAL :: tsolsrf(klon,nbsrf), qsolsrf(klon,nbsrf),snsrf(klon,nbsrf)
40	REAL :: albe(klon,nbsrf), evap(klon,nbsrf)
41	REAL :: alblw(klon,nbsrf)
42	REAL :: tsoil(klon,nsoilmx,nbsrf)
43	REAL :: radsol(klon),rain_fall(klon), snow_fall(klon)
44	REAL :: solsw(klon), sollw(klon), fder(klon)
45	cIM "slab" ocean
46	REAL :: tslab(klon), seaice(klon)
47	REAL :: frugs(klon,nbsrf), agesno(klon,nbsrf)
48	REAL :: rugmer(klon)
49	REAL :: zmea(iip1jjp1), zstd(iip1jjp1)
50	REAL :: zsig(iip1jjp1), zgam(iip1jjp1), zthe(iip1*jjp1)
51	REAL :: zpic(iip1jjp1), zval(iip1jjp1), rugsrel(iip1*jjp1)
52	REAL :: qd(iip1, jjp1, llm)
53	REAL :: pctsrf(klon, nbsrf)
54	REAL :: t_ancien(klon,klev), q_ancien(klon,klev) !
55	REAL :: run_off_lic_0(klon)
56	real :: clwcon(klon,klev),rnebcon(klon,klev),ratqs(klon,klev)
57	! declarations pour lecture glace de mer
58	REAL :: rugv(klon)
59	INTEGER :: iml_lic, jml_lic, llm_tmp, ttm_tmp, iret
60	INTEGER :: itaul(1), fid
61	REAL :: lev(1), date
62	REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_lic, lat_lic
63	REAL, ALLOCATABLE, DIMENSION(:) :: dlon_lic, dlat_lic
64	REAL, ALLOCATABLE, DIMENSION (:,:) :: fraclic
65	REAL :: flic_tmp(iip1, jjp1)
66	REAL :: champint(iim, jjp1)
67	!
68
69	CHARACTER*80 :: varname
70	!
71	INTEGER :: i,j, ig, l, ji,ii1,ii2
72	INTEGER :: nq
73	REAL :: xpi
74	!
75	REAL :: alpha(iip1,jjp1,llm),beta(iip1,jjp1,llm)
76	REAL :: pk(iip1,jjp1,llm), pls(iip1,jjp1,llm), pks(ip1jmp1)
77	REAL :: workvar(iip1,jjp1,llm)
78	!
79	REAL :: prefkap, unskap
80	!
81	real :: time_step,t_ops,t_wrt
82
83	#include "comdissnew.h"
84	#include "control.h"
85	#include "serre.h"
86	#include "clesphys.h"
87
88	INTEGER :: longcles
89	PARAMETER ( longcles = 20 )
90	REAL :: clesphy0 ( longcles )
91	REAL :: p(iip1,jjp1,llm)
92	INTEGER :: itau, iday
93	REAL :: masse(iip1,jjp1,llm)
94	REAL :: xpn,xps,xppn(iim),xpps(iim)
95	real :: time
96	REAL :: phi(ip1jmp1,llm)
97	REAL :: pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
98	REAL :: w(ip1jmp1,llm)
99	REAL ::phystep
100	INTEGER :: radpas
101	real zrel(iip1*jjp1),chmin,chmax
102
103	CHARACTER*80 :: visu_file
104	INTEGER :: visuid
105
106	! pour la lecture du fichier masque ocean
107	integer :: nid_o2a
108	logical :: couple = .false.
109	INTEGER :: iml_omask, jml_omask
110	REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_omask, lat_omask
111	REAL, ALLOCATABLE, DIMENSION(:) :: dlon_omask, dlat_omask
112	REAL, ALLOCATABLE, DIMENSION (:,:) :: ocemask, ocetmp
113	real, dimension(klon) :: ocemask_fi
114	integer :: isst(klon-2)
115	real zx_tmp_2d(iim,jjp1)
116	!
117	! Constantes
118	!
119	pi = 4. * ATAN(1.)
120	rad = 6371229.
121	omeg = 4.* ASIN(1.)/(24.*3600.)
122	g = 9.8
123	daysec = 86400.
124	kappa = 0.2857143
125	cpp = 1004.70885
126	!
127	preff = 101325.
128	unskap = 1./kappa
129	!
130	jmp1 = jjm + 1
131	!
132	! Construct a grid
133	!
134
135	! CALL defrun_new(99,.TRUE.,clesphy0)
136	CALL conf_gcm( 99, .TRUE. , clesphy0 )
137
138	dtvr = daysec/FLOAT(day_step)
139	print*,'dtvr',dtvr
140
141	CALL inicons0()
142	CALL inigeom()
143	!
144	CALL inifilr()
145	!
146	latfi(1) = ASIN(1.0)
147	DO j = 2, jjm
148	DO i = 1, iim
149	latfi((j-2)*iim+1+i)= rlatu(j)
150	ENDDO
151	ENDDO
152	latfi(klon) = - ASIN(1.0)
153	!
154	lonfi(1) = 0.0
155	DO j = 2, jjm
156	DO i = 1, iim
157	lonfi((j-2)*iim+1+i) = rlonv(i)
158	ENDDO
159	ENDDO
160	lonfi(klon) = 0.0
161	!
162	xpi = 2.0 * ASIN(1.0)
163	DO ig = 1, klon
164	latfi(ig) = latfi(ig) * 180.0 / xpi
165	lonfi(ig) = lonfi(ig) * 180.0 / xpi
166	ENDDO
167	!
168
169
170	C
171	C En cas de simulation couplee, lecture du masque ocean issu du modele ocean
172	C utilise pour calculer les poids et pour assurer l'adequation entre les
173	C fractions d'ocean vu par l'atmosphere et l'ocean. Sinon, on cree le masque
174	C a partir du fichier relief
175	C
176
177	write(,)'Essai de lecture masque ocean'
178	iret = nf_open("o2a.nc", NF_NOWRITE, nid_o2a)
179	if (iret .ne. 0) then
180	write(,)'ATTENTION!! pas de fichier o2a.nc trouve'
181	write(,)'Run force'
182	varname = 'masque'
183	masque(:,:) = 0.0
184	CALL startget(varname, iip1, jjp1, rlonv, rlatu, masque, 0.0,
185	, jjm ,rlonu,rlatv , interbar )
186	WRITE(,) 'MASQUE construit : Masque'
187	WRITE(*,'(97I1)') nINT(masque(:,:))
188	call gr_dyn_fi(1, iip1, jjp1, klon, masque, zmasq)
189	WHERE (zmasq(1 : klon) .LT. EPSFRA)
190	zmasq(1 : klon) = 0.
191	END WHERE
192	WHERE (1. - zmasq(1 : klon) .LT. EPSFRA)
193	zmasq(1 : klon) = 1.
194	END WHERE
195	else
196	couple = .true.
197	iret = nf_close(nid_o2a)
198	call flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp
199	$ , nid_o2a)
200	if (iml_omask /= iim .or. jml_omask /= jjp1) then
201	write(,)'Dimensions non compatibles pour masque ocean'
202	write(,)'iim = ',iim,' iml_omask = ',iml_omask
203	write(,)'jjp1 = ',jjp1,' jml_omask = ',jml_omask
204	stop
205	endif
206	ALLOCATE(lat_omask(iml_omask, jml_omask), stat=iret)
207	ALLOCATE(lon_omask(iml_omask, jml_omask), stat=iret)
208	ALLOCATE(dlon_omask(iml_omask), stat=iret)
209	ALLOCATE(dlat_omask(jml_omask), stat=iret)
210	ALLOCATE(ocemask(iml_omask, jml_omask), stat=iret)
211	ALLOCATE(ocetmp(iml_omask, jml_omask), stat=iret)
212	CALL flinopen("o2a.nc", .FALSE., iml_omask, jml_omask, llm_tmp
213	$ , lon_omask, lat_omask, lev, ttm_tmp, itaul, date, dt, fid)
214	CALL flinget(fid, 'OceMask', iml_omask, jml_omask, llm_tmp,
215	$ ttm_tmp, 1, 1, ocetmp)
216	CALL flinclo(fid)
217	dlon_omask(1 : iml_omask) = lon_omask(1 : iml_omask, 1)
218	dlat_omask(1 : jml_omask) = lat_omask(1 , 1 : jml_omask)
219	ocemask = ocetmp
220	if (dlat_omask(1) < dlat_omask(jml_omask)) then
221	do j = 1, jml_omask
222	ocemask(:,j) = ocetmp(:,jml_omask-j+1)
223	enddo
224	endif
225	C
226	C passage masque ocean a la grille physique
227	C
228	write(,)'ocemask '
229	write(*,'(96i1)')int(ocemask)
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	call gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque)
241
242	varname = 'relief'
243	! This line needs to be replaced by a call to restget to get the values in the restart file
244	orog(:,:) = 0.0
245	CALL startget(varname, iip1, jjp1, rlonv, rlatu, orog, 0.0 ,
246	, jjm ,rlonu,rlatv , interbar, masque )
247	!
248	WRITE(,) 'OUT OF GET VARIABLE : Relief'
249	! WRITE(*,'(49I1)') INT(orog(:,:))
250	!
251	varname = 'rugosite'
252	! This line needs to be replaced by a call to restget to get the values in the restart file
253	rugo(:,:) = 0.0
254	CALL startget(varname, iip1, jjp1, rlonv, rlatu, rugo, 0.0 ,
255	, jjm, rlonu,rlatv , interbar )
256	!
257	WRITE(,) 'OUT OF GET VARIABLE : Rugosite'
258	! WRITE(,'(49I1)') INT(rugo(:,:)10)
259	!
260	C
261	C on initialise les sous surfaces
262	C
263	pctsrf=0.
264	c
265	varname = 'psol'
266	psol(:,:) = 0.0
267	CALL startget(varname, iip1, jjp1, rlonv, rlatu, psol, 0.0 ,
268	, jjm ,rlonu,rlatv , interbar )
269	!
270	! Compute here the pressure on the intermediate levels. One would expect that this is available in the GCM
271	! anyway.
272	!
273	! WRITE(,) 'PSOL :', psol(10,20)
274	! WRITE(,) ap(:), bp(:)
275	CALL pression(ip1jmp1, ap, bp, psol, p3d)
276	! WRITE(,) 'P3D :', p3d(10,20,:)
277	CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, workvar)
278	! WRITE(,) 'PK:', pk(10,20,:)
279	!
280	!
281	!
282	prefkap = preff ** kappa
283	! WRITE(,) 'unskap, cpp, preff :', unskap, cpp, preff
284	DO l = 1, llm
285	DO j=1,jjp1
286	DO i =1, iip1
287	pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
288	ENDDO
289	ENDDO
290	ENDDO
291	!
292	! WRITE(,) 'PLS :', pls(10,20,:)
293	!
294	varname = 'surfgeo'
295	phis(:,:) = 0.0
296	CALL startget(varname, iip1, jjp1, rlonv, rlatu, phis, 0.0 ,
297	, jjm ,rlonu,rlatv, interbar )
298	!
299	varname = 'u'
300	uvent(:,:,:) = 0.0
301	CALL startget(varname, iip1, jjp1, rlonu, rlatu, llm, pls,
302	. workvar, uvent, 0.0, jjm ,rlonv, rlatv, interbar )
303	!
304	varname = 'v'
305	vvent(:,:,:) = 0.0
306	CALL startget(varname, iip1, jjm, rlonv, rlatv, llm, pls,
307	. workvar, vvent, 0.0, jjp1, rlonu, rlatu, interbar )
308	!
309	varname = 't'
310	t3d(:,:,:) = 0.0
311	CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
312	. workvar, t3d, 0.0 , jjm, rlonu, rlatv , interbar )
313	!
314	WRITE(,) 'T3D min,max:',minval(t3d(:,:,:)),
315	. maxval(t3d(:,:,:))
316	varname = 'tpot'
317	tpot(:,:,:) = 0.0
318	CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
319	. pk, tpot, 0.0 , jjm, rlonu, rlatv , interbar )
320	!
321	WRITE(,) 'T3D min,max:',minval(t3d(:,:,:)),
322	. maxval(t3d(:,:,:))
323	WRITE(,) 'PLS min,max:',minval(pls(:,:,:)),
324	. maxval(pls(:,:,:))
325
326	c Calcul de l'humidite a saturation
327	print*,'avant q_sat'
328	call q_sat(llmjjp1iip1,t3d,pls,qsat)
329	print*,'apres q_sat'
330
331	WRITE(,) 'QSAT min,max:',minval(qsat(:,:,:)),
332	. maxval(qsat(:,:,:))
333	!
334	WRITE(,) 'QSAT :', qsat(10,20,:)
335	!
336	varname = 'q'
337	qd(:,:,:) = 0.0
338	q3d(:,:,:,:) = 0.0
339	WRITE(,) 'QSAT min,max:',minval(qsat(:,:,:)),
340	. maxval(qsat(:,:,:))
341	CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
342	. qsat, qd, 0.0, jjm, rlonu, rlatv , interbar )
343	q3d(:,:,:,1) = qd(:,:,:)
344	!
345	varname = 'tsol'
346	! This line needs to be replaced by a call to restget to get the values in the restart file
347	tsol(:) = 0.0
348	CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, tsol, 0.0,
349	. jjm, rlonu, rlatv , interbar )
350	!
351	WRITE(,) 'TSOL construit :'
352	! WRITE(*,'(48I3)') INT(TSOL(2:klon)-273)
353	!
354	varname = 'qsol'
355	qsol(:) = 0.0
356	CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, qsol, 0.0,
357	. jjm, rlonu, rlatv , interbar )
358	!
359	varname = 'snow'
360	sn(:) = 0.0
361	CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, sn, 0.0,
362	. jjm, rlonu, rlatv , interbar )
363	!
364	varname = 'rads'
365	radsol(:) = 0.0
366	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,radsol,0.0,
367	. jjm, rlonu, rlatv , interbar )
368	!
369	cIM "slab" ocean
370	varname = 'tslab'
371	tslab(:) = 0.0
372	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,tslab,0.0,
373	. jjm, rlonu, rlatv , interbar )
374	c
375	varname = 'seaice'
376	seaice(:) = 0.0
377	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,seaice,0.0,
378	. jjm, rlonu, rlatv , interbar )
379	!
380	varname = 'rugmer'
381	rugmer(:) = 0.0
382	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,rugmer,0.0,
383	. jjm, rlonu, rlatv , interbar )
384	!
385	! varname = 'agesno'
386	! agesno(:) = 0.0
387	! CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,agesno,0.0,
388	! . jjm, rlonu, rlatv , interbar )
389
390	varname = 'zmea'
391	zmea(:) = 0.0
392	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zmea,0.0,
393	. jjm, rlonu, rlatv , interbar )
394
395	varname = 'zstd'
396	zstd(:) = 0.0
397	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zstd,0.0,
398	. jjm, rlonu, rlatv , interbar )
399	varname = 'zsig'
400	zsig(:) = 0.0
401	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zsig,0.0,
402	. jjm, rlonu, rlatv , interbar )
403	varname = 'zgam'
404	zgam(:) = 0.0
405	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zgam,0.0,
406	. jjm, rlonu, rlatv , interbar )
407	varname = 'zthe'
408	zthe(:) = 0.0
409	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zthe,0.0,
410	. jjm, rlonu, rlatv , interbar )
411	varname = 'zpic'
412	zpic(:) = 0.0
413	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zpic,0.0,
414	. jjm, rlonu, rlatv , interbar )
415	varname = 'zval'
416	zval(:) = 0.0
417	CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zval,0.0,
418	. jjm, rlonu, rlatv , interbar )
419	c
420	rugsrel(:) = 0.0
421	IF(ok_orodr) THEN
422	DO i = 1, iip1* jjp1
423	rugsrel(i) = MAX( 1.e-05, zstd(i)* zsig(i) /2. )
424	ENDDO
425	ENDIF
426
427
428	C
429	C lecture du fichier glace de terre pour fixer la fraction de terre
430	C et de glace de terre
431	C
432	CALL flininfo("landiceref.nc", iml_lic, jml_lic,llm_tmp, ttm_tmp
433	$ , fid)
434	ALLOCATE(lat_lic(iml_lic, jml_lic), stat=iret)
435	ALLOCATE(lon_lic(iml_lic, jml_lic), stat=iret)
436	ALLOCATE(dlon_lic(iml_lic), stat=iret)
437	ALLOCATE(dlat_lic(jml_lic), stat=iret)
438	ALLOCATE(fraclic(iml_lic, jml_lic), stat=iret)
439	CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp
440	$ , lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid)
441	CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp
442	$ , 1, 1, fraclic)
443	CALL flinclo(fid)
444	C
445	C interpolation sur la grille T du modele
446	C
447	WRITE(,) 'dimensions de landice iml_lic, jml_lic : ',
448	$ iml_lic, jml_lic
449	c
450	C sil les coordonnees sont en degres, on les transforme
451	C
452	IF( MAXVAL( lon_lic(:,:) ) .GT. 2.0 * asin(1.0) ) THEN
453	lon_lic(:,:) = lon_lic(:,:) * 2.0* ASIN(1.0) / 180.
454	ENDIF
455	IF( maxval( lat_lic(:,:) ) .GT. 2.0 * asin(1.0)) THEN
456	lat_lic(:,:) = lat_lic(:,:) * 2.0 * asin(1.0) / 180.
457	ENDIF
458
459	dlon_lic(1 : iml_lic) = lon_lic(1 : iml_lic, 1)
460	dlat_lic(1 : jml_lic) = lat_lic(1 , 1 : jml_lic)
461	C
462	CALL grille_m(iml_lic, jml_lic, dlon_lic, dlat_lic, fraclic
463	$ ,iim, jjp1,
464	$ rlonv, rlatu, flic_tmp(1 : iim, 1 : jjp1))
465	c$$$ flic_tmp(1 : iim, 1 : jjp1) = champint(1: iim, 1 : jjp1)
466	flic_tmp(iip1, 1 : jjp1) = flic_tmp(1 , 1 : jjp1)
467	C
468	C passage sur la grille physique
469	C
470	CALL gr_dyn_fi(1, iip1, jjp1, klon, flic_tmp,
471	$ pctsrf(1:klon, is_lic))
472	C adequation avec le maque terre/mer
473	c zmasq(157) = 0.
474	WHERE (pctsrf(1 : klon, is_lic) .LT. EPSFRA )
475	pctsrf(1 : klon, is_lic) = 0.
476	END WHERE
477	WHERE (zmasq( 1 : klon) .LT. EPSFRA)
478	pctsrf(1 : klon, is_lic) = 0.
479	END WHERE
480	pctsrf(1 : klon, is_ter) = zmasq(1 : klon)
481	DO ji = 1, klon
482	IF (zmasq(ji) .GT. EPSFRA) THEN
483	IF ( pctsrf(ji, is_lic) .GE. zmasq(ji)) THEN
484	pctsrf(ji, is_lic) = zmasq(ji)
485	pctsrf(ji, is_ter) = 0.
486	ELSE
487	pctsrf(ji,is_ter) = zmasq(ji) - pctsrf(ji, is_lic)
488	IF (pctsrf(ji,is_ter) .LT. EPSFRA) THEN
489	pctsrf(ji,is_ter) = 0.
490	pctsrf(ji, is_lic) = zmasq(ji)
491	ENDIF
492	ENDIF
493	ENDIF
494	END DO
495	C
496	C sous surface ocean et glace de mer (pour demarrer on met glace de mer a 0)
497	C
498	pctsrf(1 : klon, is_oce) = (1. - zmasq(1 : klon))
499
500
501	WHERE (pctsrf(1 : klon, is_oce) .LT. EPSFRA)
502	pctsrf(1 : klon, is_oce) = 0.
503	END WHERE
504
505	if (couple) pctsrf(1 : klon, is_oce) = ocemask_fi(1 : klon)
506
507	isst = 0
508	where (pctsrf(2:klon-1,is_oce) >0.) isst = 1
509	C
510	C verif que somme des sous surface = 1
511	C
512	ji=count(abs(sum(pctsrf(1 : klon, 1 : nbsrf), dim = 2) - 1.0 )
513	$ .GT. EPSFRA)
514	IF (ji .NE. 0) THEN
515	WRITE(,) 'pb repartition sous maille pour ',ji,' points'
516	ENDIF
517
518	! where (pctsrf(1:klon, is_ter) >= .5)
519	! pctsrf(1:klon, is_ter) = 1.
520	! pctsrf(1:klon, is_oce) = 0.
521	! pctsrf(1:klon, is_sic) = 0.
522	! pctsrf(1:klon, is_lic) = 0.
523	! zmasq = 1.
524	! endwhere
525	! where (pctsrf(1:klon, is_lic) >= .5)
526	! pctsrf(1:klon, is_ter) = 0.
527	! pctsrf(1:klon, is_oce) = 0.
528	! pctsrf(1:klon, is_sic) = 0.
529	! pctsrf(1:klon, is_lic) = 1.
530	! zmasq = 1.
531	! endwhere
532	! where (pctsrf(1:klon, is_oce) >= .5)
533	! pctsrf(1:klon, is_ter) = 0.
534	! pctsrf(1:klon, is_oce) = 1.
535	! pctsrf(1:klon, is_sic) = 0.
536	! pctsrf(1:klon, is_lic) = 0.
537	! zmasq = 0.
538	! endwhere
539	! where (pctsrf(1:klon, is_sic) >= .5)
540	! pctsrf(1:klon, is_ter) = 0.
541	! pctsrf(1:klon, is_oce) = 0.
542	! pctsrf(1:klon, is_sic) = 1.
543	! pctsrf(1:klon, is_lic) = 0.
544	! zmasq = 0.
545	! endwhere
546	! call gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque)
547	C
548	C verif que somme des sous surface = 1
549	C
550	! ji=count( (abs( sum(pctsrf(1 : klon, 1 : nbsrf), dim = 2)) - 1.0 )
551	! $ .GT. EPSFRA)
552	! IF (ji .NE. 0) THEN
553	! WRITE(,) 'pb repartition sous maille pour ',ji,' points'
554	! ENDIF
555
556	CALL gr_fi_ecrit(1,klon,iim,jjp1,zmasq,zx_tmp_2d)
557	write(,)'zmasq = '
558	write(*,'(96i1)')nint(zx_tmp_2d)
559	call gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque)
560	WRITE(,) 'MASQUE construit : Masque'
561	WRITE(*,'(97I1)') nINT(masque(:,:))
562
563
564
565	C Calcul intermediaire
566	c
567	CALL massdair( p3d, masse )
568	c
569
570	print *,' ALPHAX ',alphax
571
572	DO l = 1, llm
573	DO i = 1, iim
574	xppn(i) = aire( i, 1 ) * masse( i , 1 , l )
575	xpps(i) = aire( i,jjp1 ) * masse( i , jjp1 , l )
576	ENDDO
577	xpn = SUM(xppn)/apoln
578	xps = SUM(xpps)/apols
579	DO i = 1, iip1
580	masse( i , 1 , l ) = xpn
581	masse( i , jjp1 , l ) = xps
582	ENDDO
583	ENDDO
584	q3d(iip1,:,:,:) = q3d(1,:,:,:)
585	phis(iip1,:) = phis(1,:)
586
587	C init pour traceurs
588	call iniadvtrac(nq)
589	C Ecriture
590	CALL inidissip( lstardis, nitergdiv, nitergrot, niterh ,
591	* tetagdiv, tetagrot , tetatemp )
592	print*,'sortie inidissip'
593	itau = 0
594	itau_dyn = 0
595	itau_phy = 0
596	iday = dayref +itau/day_step
597	time = FLOAT(itau-(iday-dayref)*day_step)/day_step
598	c
599	IF(time.GT.1) THEN
600	time = time - 1
601	iday = iday + 1
602	ENDIF
603	day_ref = dayref
604	annee_ref = anneeref
605
606	CALL geopot ( ip1jmp1, tpot , pk , pks, phis , phi )
607	print*,'sortie geopot'
608
609	CALL caldyn0 ( itau,uvent,vvent,tpot,psol,masse,pk,phis ,
610	* phi,w, pbaru,pbarv,time+iday-dayref )
611	print*,'sortie caldyn0'
612	CALL dynredem0("start.nc",dayref,phis,nqmx)
613	print*,'sortie dynredem0'
614	CALL dynredem1("start.nc",0.0,vvent,uvent,tpot,q3d,nqmx,masse ,
615	. psol)
616	print*,'sortie dynredem1'
617	C
618	C Ecriture etat initial physique
619	C
620	write(,)'phystep ',dtvr,iphysiq,nbapp_rad
621	phystep = dtvr * FLOAT(iphysiq)
622	radpas = NINT (86400./phystep/ FLOAT(nbapp_rad) )
623	write(,)'phystep =', phystep, radpas
624	cIM : lecture de co2_ppm & solaire ds physiq.def
625	c co2_ppm = 348.0
626	c solaire = 1365.0
627
628	c
629	c Initialisation
630	c tsol, qsol, sn,albe, evap,tsoil,rain_fall, snow_fall,solsw, sollw,frugs
631	c
632	tsolsrf(:,is_ter) = tsol
633	tsolsrf(:,is_lic) = tsol
634	tsolsrf(:,is_oce) = tsol
635	tsolsrf(:,is_sic) = tsol
636	snsrf(:,is_ter) = sn
637	snsrf(:,is_lic) = sn
638	snsrf(:,is_oce) = sn
639	snsrf(:,is_sic) = sn
640	albe(:,is_ter) = 0.08
641	albe(:,is_lic) = 0.6
642	albe(:,is_oce) = 0.5
643	albe(:,is_sic) = 0.6
644	alblw = albe
645	evap(:,:) = 0.
646	qsolsrf(:,is_ter) = 150
647	qsolsrf(:,is_lic) = 150
648	qsolsrf(:,is_oce) = 150.
649	qsolsrf(:,is_sic) = 150.
650	do i = 1, nbsrf
651	do j = 1, nsoilmx
652	tsoil(:,j,i) = tsol
653	enddo
654	enddo
655	rain_fall = 0.; snow_fall = 0.
656	solsw = 165.
657	sollw = -53.
658	t_ancien = 273.15
659	q_ancien = 0.
660	agesno = 0.
661	cIM "slab" ocean
662	tslab(1:klon) = tsolsrf(1:klon,is_oce)
663	seaice = 0. !kg m-2
664	c
665	frugs(1:klon,is_oce) = rugmer(1:klon)
666	frugs(1:klon,is_ter) = MAX(1.0e-05, zstd(1:klon)*zsig(1:klon)/2.0)
667	frugs(1:klon,is_lic) = MAX(1.0e-05, zstd(1:klon)*zsig(1:klon)/2.0)
668	frugs(1:klon,is_sic) = 0.001
669	fder = 0.0
670	clwcon = 0.0
671	rnebcon = 0.0
672	ratqs = 0.0
673	run_off_lic_0 = 0.0
674
675	cIM call phyredem("startphy.nc",phystep,radpas, co2_ppm, solaire,
676	call phyredem("startphy.nc",phystep,radpas,
677	$ latfi, lonfi, pctsrf, tsolsrf, tsoil, tslab, seaice,
678	$ qsolsrf, qsol,
679	$ snsrf,
680	$ albe, alblw, evap, rain_fall, snow_fall, solsw, sollw, fder,
681	$ radsol, frugs, agesno,
682	$ zmea, zstd, zsig, zgam, zthe, zpic, zval, rugsrel,
683	$ t_ancien, q_ancien, rnebcon, ratqs, clwcon,
684	$ run_off_lic_0)
685	print*,'sortie phyredem'
686
687	C Sortie Visu pour les champs dynamiques
688	if (1.eq.0 ) then
689	print*,'sortie visu'
690	time_step = 1.
691	t_ops = 2.
692	t_wrt = 2.
693	itau = 2.
694	visu_file='Etat0_visu.nc'
695	CALL initdynav(visu_file,dayref,anneeref,time_step,
696	. t_ops, t_wrt, nqmx, visuid)
697	CALL writedynav(visuid, nqmx, itau,vvent ,
698	. uvent,tpot,pk,phi,q3d,masse,psol,phis)
699	else
700	print*,'CCCCCCCCCCCCCCCCCC REACTIVER SORTIE VISU DANS ETAT0'
701	endif
702	print*,'entree histclo'
703	CALL histclo
704	RETURN
705	!
706	END SUBROUTINE etat0_netcdf

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

Download in other formats:

Original Format