Context Navigation

etat0_netcdf.F @ 1174

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

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

Download in other formats:

Original Format