Context Navigation

etat0_netcdf.F @ 1198

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

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

Download in other formats:

Original Format