Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

etat0_netcdf.F @ 371

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

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

Download in other formats:

Original Format