Context Navigation

etat0_netcdf.F @ 1256

Last change on this file since 1256 was 1227, checked in by jghattas, 15 years ago

Inclusion d'un premier version du cycle de carbon dans LMDZ. Attention

!! Il s'agit d'un version ou les nouveaux cles cycle_carbon_tr et
cycle_carbon_cpl ne sont pas teste. Avec les ancinenes parametres le
modele donne les memes resultats qu'avant. L'interface avec ORCHIDEE n'a
pas encore etait modifie.

physiq.F, phys_cal_mod.F90 : ajout d'un nouveau module qui contient qq parametres pour le calendrier et le pas de temps acutelle de la physiq. Ce module pourrait etre elargie plus tard / LF + JG

infotrac.F90 : les noms du traceurs peut prendre un nom plus long (15 caracteres) dans traceur.def

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

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

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

Download in other formats:

Original Format