Context Navigation

etat0_netcdf.F @ 971

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

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

Download in other formats:

Original Format