Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

etat0_netcdf.F @ 402

Last change on this file since 402 was 402, checked in by (none), 22 years ago
This commit was manufactured by cvs2svn to create tag 'IPSL-CM4_v1'.
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 21.5 KB

Line
1	c
2	c $Header$
3	c
4	SUBROUTINE etat0_netcdf (interbar, masque, pctsrf)
5
6	USE startvar
7	USE ioipsl
8	!
9	IMPLICIT NONE
10	!
11	#include "netcdf.inc"
12	#include "dimensions.h"
13	#include "paramet.h"
14	!
15	!
16	! INTEGER, PARAMETER :: KIDIA=1, KFDIA=iim*(jjm-1)+2,
17	! .KLON=KFDIA-KIDIA+1,KLEV=llm
18	!
19	#include "comgeom2.h"
20	#include "comvert.h"
21	#include "comconst.h"
22	#include "indicesol.h"
23	#include "dimphy.h"
24	#include "dimsoil.h"
25	#include "temps.h"
26	!
27	LOGICAL interbar
28	REAL :: latfi(klon), lonfi(klon)
29	REAL :: orog(iip1,jjp1), rugo(iip1,jjp1), masque(iip1,jjp1),
30	. psol(iip1, jjp1), phis(iip1, jjp1)
31	REAL :: p3d(iip1, jjp1, llm+1)
32	REAL :: uvent(iip1, jjp1, llm)
33	REAL :: vvent(iip1, jjm, llm)
34	REAL :: t3d(iip1, jjp1, llm), tpot(iip1, jjp1, llm)
35	REAL :: q3d(iip1, jjp1, llm,nqmx), qsat(iip1, jjp1, llm)
36	REAL :: tsol(klon), qsol(klon), sn(klon)
37	REAL :: tsolsrf(klon,nbsrf), qsolsrf(klon,nbsrf),snsrf(klon,nbsrf)
38	REAL :: albe(klon,nbsrf), evap(klon,nbsrf)
39	REAL :: tsoil(klon,nsoilmx,nbsrf)
40	REAL :: radsol(klon),rain_fall(klon), snow_fall(klon)
41	REAL :: solsw(klon), sollw(klon), fder(klon)
42	REAL :: deltat(klon), frugs(klon,nbsrf), agesno(klon,nbsrf)
43	REAL :: rugmer(klon)
44	REAL :: zmea(iip1jjp1), zstd(iip1jjp1)
45	REAL :: zsig(iip1jjp1), zgam(iip1jjp1), zthe(iip1*jjp1)
46	REAL :: zpic(iip1jjp1), zval(iip1jjp1), rugsrel(iip1*jjp1)
47	REAL :: qd(iip1, jjp1, llm)
48	REAL :: pctsrf(klon, nbsrf)
49	REAL :: t_ancien(klon,klev), q_ancien(klon,klev) !
50	real :: clwcon(klon,klev),rnebcon(klon,klev),ratqs(klon,klev)
51	! declarations pour lecture glace de mer
52	REAL :: rugv(klon)
53	INTEGER :: iml_lic, jml_lic, llm_tmp, ttm_tmp, iret
54	INTEGER :: itaul(1), fid
55	REAL :: lev(1), date
56	REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_lic, lat_lic
57	REAL, ALLOCATABLE, DIMENSION(:) :: dlon_lic, dlat_lic
58	REAL, ALLOCATABLE, DIMENSION (:,:) :: fraclic
59	REAL :: flic_tmp(iip1, jjp1)
60	REAL :: champint(iim, jjp1)
61	!
62
63	CHARACTER*80 :: varname
64	!
65	INTEGER :: i,j, ig, l, ji,ii1,ii2
66	REAL :: xpi
67	!
68	REAL :: alpha(iip1,jjp1,llm),beta(iip1,jjp1,llm)
69	REAL :: pk(iip1,jjp1,llm), pls(iip1,jjp1,llm), pks(ip1jmp1)
70	REAL :: workvar(iip1,jjp1,llm)
71	!
72	REAL :: prefkap, unskap
73	!
74	REAL :: q_sat
75	EXTERNAL q_sat
76	real :: time_step,t_ops,t_wrt
77
78	#include "comdissnew.h"
79	#include "control.h"
80	#include "serre.h"
81	#include "clesphys.h"
82
83	INTEGER :: longcles
84	PARAMETER ( longcles = 20 )
85	REAL :: clesphy0 ( longcles )
86	REAL :: p(iip1,jjp1,llm)
87	INTEGER :: itau, iday
88	REAL :: masse(iip1,jjp1,llm)
89	REAL :: xpn,xps,xppn(iim),xpps(iim)
90	real :: time
91	REAL :: phi(ip1jmp1,llm)
92	REAL :: pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
93	REAL :: w(ip1jmp1,llm)
94	REAL ::phystep,co2_ppm,solaire
95	INTEGER :: radpas
96	real zrel(iip1*jjp1),chmin,chmax
97
98	CHARACTER*80 :: visu_file
99	INTEGER :: visuid
100
101	! pour la lecture du fichier masque ocean
102	integer :: nid_o2a
103	logical :: couple = .false.
104	INTEGER :: iml_omask, jml_omask
105	REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_omask, lat_omask
106	REAL, ALLOCATABLE, DIMENSION(:) :: dlon_omask, dlat_omask
107	REAL, ALLOCATABLE, DIMENSION (:,:) :: ocemask, ocetmp
108	real, dimension(klon) :: ocemask_fi
109	integer :: isst(klon-2)
110	real zx_tmp_2d(iim,jjp1)
111	!
112	! Constantes
113	!
114	pi = 4. * ATAN(1.)
115	rad = 6371229.
116	omeg = 4.* ASIN(1.)/(24.*3600.)
117	g = 9.8
118	daysec = 86400.
119	kappa = 0.2857143
120	cpp = 1004.70885
121	!
122	preff = 101325.
123	unskap = 1./kappa
124	!
125	jmp1 = jjm + 1
126	!
127	! Construct a grid
128	!
129
130	! CALL defrun_new(99,.TRUE.,clesphy0)
131	CALL conf_gcm( 99, .TRUE. , clesphy0 )
132
133	dtvr = daysec/FLOAT(day_step)
134	print*,'dtvr',dtvr
135
136	CALL inicons0()
137	CALL inigeom()
138	!
139	CALL inifilr()
140	!
141	latfi(1) = ASIN(1.0)
142	DO j = 2, jjm
143	DO i = 1, iim
144	latfi((j-2)*iim+1+i)= rlatu(j)
145	ENDDO
146	ENDDO
147	latfi(klon) = - ASIN(1.0)
148	!
149	lonfi(1) = 0.0
150	DO j = 2, jjm
151	DO i = 1, iim
152	lonfi((j-2)*iim+1+i) = rlonv(i)
153	ENDDO
154	ENDDO
155	lonfi(klon) = 0.0
156	!
157	xpi = 2.0 * ASIN(1.0)
158	DO ig = 1, klon
159	latfi(ig) = latfi(ig) * 180.0 / xpi
160	lonfi(ig) = lonfi(ig) * 180.0 / xpi
161	ENDDO
162	!
163
164
165	C
166	C En cas de simulation couplee, lecture du masque ocean issu du modele ocean
167	C utilise pour calculer les poids et pour assurer l'adequation entre les
168	C fractions d'ocean vu par l'atmosphere et l'ocean. Sinon, on cree le masque
169	C a partir du fichier relief
170	C
171
172	write(,)'Essai de lecture masque ocean'
173	iret = nf_open("o2a.nc", NF_NOWRITE, nid_o2a)
174	if (iret .ne. 0) then
175	write(,)'ATTENTION!! pas de fichier o2a.nc trouve'
176	write(,)'Run force'
177	varname = 'masque'
178	masque(:,:) = 0.0
179	CALL startget(varname, iip1, jjp1, rlonv, rlatu, masque, 0.0,
180	, jjm ,rlonu,rlatv , interbar )
181	WRITE(,) 'MASQUE construit : Masque'
182	WRITE(*,'(97I1)') nINT(masque(:,:))
183	call gr_dyn_fi(1, iip1, jjp1, klon, masque, zmasq)
184	WHERE (zmasq(1 : klon) .LT. EPSFRA)
185	zmasq(1 : klon) = 0.
186	END WHERE
187	WHERE (1. - zmasq(1 : klon) .LT. EPSFRA)
188	zmasq(1 : klon) = 1.
189	END WHERE
190	else
191	couple = .true.
192	iret = nf_close(nid_o2a)
193	call flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp
194	$ , nid_o2a)
195	if (iml_omask /= iim .or. jml_omask /= jjp1) then
196	write(,)'Dimensions non compatibles pour masque ocean'
197	write(,)'iim = ',iim,' iml_omask = ',iml_omask
198	write(,)'jjp1 = ',jjp1,' jml_omask = ',jml_omask
199	stop
200	endif
201	ALLOCATE(lat_omask(iml_omask, jml_omask), stat=iret)
202	ALLOCATE(lon_omask(iml_omask, jml_omask), stat=iret)
203	ALLOCATE(dlon_omask(iml_omask), stat=iret)
204	ALLOCATE(dlat_omask(jml_omask), stat=iret)
205	ALLOCATE(ocemask(iml_omask, jml_omask), stat=iret)
206	ALLOCATE(ocetmp(iml_omask, jml_omask), stat=iret)
207	CALL flinopen("o2a.nc", .FALSE., iml_omask, jml_omask, llm_tmp
208	$ , lon_omask, lat_omask, lev, ttm_tmp, itaul, date, dt, fid)
209	CALL flinget(fid, 'OceMask', iml_omask, jml_omask, llm_tmp,
210	$ ttm_tmp, 1, 1, ocetmp)
211	CALL flinclo(fid)
212	dlon_omask(1 : iml_omask) = lon_omask(1 : iml_omask, 1)
213	dlat_omask(1 : jml_omask) = lat_omask(1 , 1 : jml_omask)
214	ocemask = ocetmp
215	if (dlat_omask(1) < dlat_omask(jml_omask)) then
216	do j = 1, jml_omask
217	ocemask(:,j) = ocetmp(:,jml_omask-j+1)
218	enddo
219	endif
220	C
221	C passage masque ocean a la grille physique
222	C
223	write(,)'ocemask '
224	write(*,'(96i1)')int(ocemask)
225	ocemask_fi(1) = ocemask(1,1)
226	do j = 2, jjm
227	do i = 1, iim
228	ocemask_fi((j-2)*iim + i + 1) = ocemask(i,j)
229	enddo
230	enddo
231	ocemask_fi(klon) = ocemask(1,jjp1)
232	zmasq = 1. - ocemask_fi
233	endif
234
235	call gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque)
236
237	varname = 'relief'
238	! This line needs to be replaced by a call to restget to get the values in the restart file
239	orog(:,:) = 0.0
240	CALL startget(varname, iip1, jjp1, rlonv, rlatu, orog, 0.0 ,
241	, jjm ,rlonu,rlatv , interbar, masque )
242	!
243	WRITE(,) 'OUT OF GET VARIABLE : Relief'
244	! WRITE(*,'(49I1)') INT(orog(:,:))
245	!
246	varname = 'rugosite'
247	! This line needs to be replaced by a call to restget to get the values in the restart file
248	rugo(:,:) = 0.0
249	CALL startget(varname, iip1, jjp1, rlonv, rlatu, rugo, 0.0 ,
250	, jjm, rlonu,rlatv , interbar )
251	!
252	WRITE(,) 'OUT OF GET VARIABLE : Rugosite'
253	! WRITE(,'(49I1)') INT(rugo(:,:)10)
254	!
255	C
256	C on initialise les sous surfaces
257	C
258	pctsrf=0.
259	c
260	varname = 'psol'
261	psol(:,:) = 0.0
262	CALL startget(varname, iip1, jjp1, rlonv, rlatu, psol, 0.0 ,
263	, jjm ,rlonu,rlatv , interbar )
264	!
265	! Compute here the pressure on the intermediate levels. One would expect that this is available in the GCM
266	! anyway.
267	!
268	! WRITE(,) 'PSOL :', psol(10,20)
269	! WRITE(,) ap(:), bp(:)
270	CALL pression(ip1jmp1, ap, bp, psol, p3d)
271	! WRITE(,) 'P3D :', p3d(10,20,:)
272	CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, workvar)
273	! WRITE(,) 'PK:', pk(10,20,:)
274	!
275	!
276	!
277	prefkap = preff ** kappa
278	! WRITE(,) 'unskap, cpp, preff :', unskap, cpp, preff
279	DO l = 1, llm
280	DO j=1,jjp1
281	DO i =1, iip1
282	pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
283	ENDDO
284	ENDDO
285	ENDDO
286	!
287	! WRITE(,) 'PLS :', pls(10,20,:)
288	!
289	varname = 'surfgeo'
290	phis(:,:) = 0.0
291	CALL startget(varname, iip1, jjp1, rlonv, rlatu, phis, 0.0 ,
292	, jjm ,rlonu,rlatv, interbar )
293	!
294	varname = 'u'
295	uvent(:,:,:) = 0.0
296	CALL startget(varname, iip1, jjp1, rlonu, rlatu, llm, pls,
297	. workvar, uvent, 0.0, jjm ,rlonv, rlatv, interbar )
298	!
299	varname = 'v'
300	vvent(:,:,:) = 0.0
301	CALL startget(varname, iip1, jjm, rlonv, rlatv, llm, pls,
302	. workvar, vvent, 0.0, jjp1, rlonu, rlatu, interbar )
303	!
304	varname = 't'
305	t3d(:,:,:) = 0.0
306	CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
307	. workvar, t3d, 0.0 , jjm, rlonu, rlatv , interbar )
308	!
309	WRITE(,) 'T3D min,max:',minval(t3d(:,:,:)),
310	. maxval(t3d(:,:,:))
311	varname = 'tpot'
312	tpot(:,:,:) = 0.0
313	CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
314	. pk, tpot, 0.0 , jjm, rlonu, rlatv , interbar )
315	!
316	WRITE(,) 'T3D min,max:',minval(t3d(:,:,:)),
317	. maxval(t3d(:,:,:))
318	WRITE(,) 'PLS min,max:',minval(pls(:,:,:)),
319	. maxval(pls(:,:,:))
320	DO l = 1, llm
321	DO j=1,jjp1
322	DO i =1, iip1-1
323	qsat(i,j,l) = q_sat(t3d(i,j,l),pls(i,j,l)/100. )
324	ENDDO
325	qsat(iip1,j,l) = qsat(1,j,l)
326	ENDDO
327	ENDDO
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 Ecriture
579
580
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	co2_ppm = 330.0
616	solaire = 1370.0
617
618	c
619	c Initialisation
620	c tsol, qsol, sn,albe, evap,tsoil,rain_fall, snow_fall,solsw, sollw,frugs
621	c
622	tsolsrf(:,is_ter) = tsol
623	tsolsrf(:,is_lic) = tsol
624	tsolsrf(:,is_oce) = tsol
625	tsolsrf(:,is_sic) = tsol
626	snsrf(:,is_ter) = sn
627	snsrf(:,is_lic) = sn
628	snsrf(:,is_oce) = sn
629	snsrf(:,is_sic) = sn
630	albe(:,is_ter) = 0.08
631	albe(:,is_lic) = 0.6
632	albe(:,is_oce) = 0.5
633	albe(:,is_sic) = 0.6
634	evap(:,:) = 0.
635	qsolsrf(:,is_ter) = 150
636	qsolsrf(:,is_lic) = 150
637	qsolsrf(:,is_oce) = 150.
638	qsolsrf(:,is_sic) = 150.
639	do i = 1, nbsrf
640	do j = 1, nsoilmx
641	tsoil(:,j,i) = tsol
642	enddo
643	enddo
644	rain_fall = 0.; snow_fall = 0.
645	solsw = 165.
646	sollw = -53.
647	t_ancien = 273.15
648	q_ancien = 0.
649	agesno = 0.
650	deltat = 0.
651	frugs(1:klon,is_oce) = rugmer(1:klon)
652	frugs(1:klon,is_ter) = MAX(1.0e-05, zstd(1:klon)*zsig(1:klon)/2.0)
653	frugs(1:klon,is_lic) = MAX(1.0e-05, zstd(1:klon)*zsig(1:klon)/2.0)
654	frugs(1:klon,is_sic) = 0.001
655	fder = 0.0
656	clwcon = 0.0
657	rnebcon = 0.0
658	ratqs = 0.0
659
660	call phyredem("startphy.nc",phystep,radpas, co2_ppm, solaire,
661	$ latfi, lonfi, pctsrf, tsolsrf, tsoil, deltat, qsolsrf, snsrf,
662	$ albe, evap, rain_fall, snow_fall, solsw, sollw, fder,
663	$ radsol, frugs, agesno,
664	$ zmea, zstd, zsig, zgam, zthe, zpic, zval, rugsrel,
665	$ t_ancien, q_ancien, rnebcon, ratqs, clwcon)
666
667	C Sortie Visu pour les champs dynamiques
668	print*,'sortie visu'
669	time_step = 1.
670	t_ops = 2.
671	t_wrt = 2.
672	itau = 2.
673	visu_file='Etat0_visu.nc'
674	CALL initdynav(visu_file,dayref,anneeref,time_step,
675	. t_ops, t_wrt, nqmx, visuid)
676	CALL writedynav(visuid, nqmx, itau,vvent ,
677	. uvent,tpot,pk,phi,q3d,masse,psol,phis)
678	print*,'entree histclo'
679	CALL histclo
680	RETURN
681	!
682	END SUBROUTINE etat0_netcdf

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

Download in other formats:

Original Format