Context Navigation

etat0_netcdf.F @ 1167

Last change on this file since 1167 was 1146, checked in by Laurent Fairhead, 16 years ago
Réintegration dans le tronc des modifications issues de la branche LMDZ-dev comprises entre la révision 1074 et 1145 Validation: une simulation de 1 jour en séquentiel sur PC donne les mêmes résultats entre la trunk et la dev LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 23.6 KB

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

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

Download in other formats:

Original Format