Context Navigation

etat0_netcdf.F @ 1056

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

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

Download in other formats:

Original Format