Context Navigation

etat0_netcdf.F @ 663

Last change on this file since 663 was 644, checked in by Laurent Fairhead, 20 years ago
Synchronisation avec tous les diagnostiques de Ionela IM Inclusion du slab ocean IM 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