Context Navigation

etat0_netcdf.F @ 1012

Last change on this file since 1012 was 1012, checked in by lsce, 16 years ago

Error in another argument list.
Added writing in double precision in case of NC_DOUBLE

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 23.2 KB

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

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

Download in other formats:

Original Format