Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

etat0_netcdf.F

Last change on this file was 1379, checked in by lguez, 15 years ago

Added optional ozone tracer with chemistry parameterized by Daniel
Cariolle. This tracer is passive: it has no influence on the rest of
the simulation.

Added variable "zmasse" in file "histrac.nc". Corrected long name of
variable "pplay" in "histrac.nc". Changed name of variable "t" to "T"
in "histrac.nc", corrected long name and unit.

In "phytrac", moved definition of "zmasse" toward the beginning of the
procedure, so that "zmasse" can be given as argument to
"traclmdz". Also added arguments "julien", "gmtime" and "xlon" to
"traclmdz". The four additional arguments are required for the ozone
tracer.

In module "traclmdz_mod", factorized declaration "implicit none" that
was in each procedure. There is now an equivalent single declaration
at the module level.

In procedure "traclmdz", removed variable "delp". Use "zmasse * rg"
instead since we now have "zmasse" as an argument.

Tests. Compilations on Brodie only, with optimization options "-debug"
and "-dev", parallelization options "none", "mpi", "omp" and
"mpi_omp", this revision and revision 1373. Run cases with and without
ozone tracer, 1 and 2 MPI processes, 1 and 2 OpenMP
threads. Comparisons of all cases are ok, except for strange
variations in variables "d_tr_cl_RN" and "d_tr_cl_PB" of file
"histrac.nc", variables "RN" and "PB" of "restart.nc", variables
"trs_RN" and "trs_PB" of "restartphy.nc". Relative variations of these
variables between cases are of order 1e-7 or less, after one day of
simulation.

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

File size: 25.4 KB

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

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

Download in other formats:

Original Format