Context Navigation

etat0_netcdf.F @ 1322

Last change on this file since 1322 was 1322, checked in by Laurent Fairhead, 15 years ago

Improvements concerning wake parametrisation (from JYG, NR, IT, with more to come).
Alp_offset is read in form physiq.def file

Améliorations à la paramétrisation des poches froides (de JYG, NR, IT, d'autres
sont à venir)
Alp_offset est rajouté à la liste des paramètres lus dans physiq.def

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

File size: 24.7 KB

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

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

Download in other formats:

Original Format