Context Navigation

etat0_netcdf.F @ 1227

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

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

Download in other formats:

Original Format