Context Navigation

etat0_netcdf.F @ 1076

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

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

Download in other formats:

Original Format