Context Navigation

etat0_netcdf.F @ 1630

Last change on this file since 1630 was 1630, checked in by Laurent Fairhead, 12 years ago

Importation initiale du rÃ©pertoire dyn3dmem

Initial import of dyn3dmem directory

File size: 25.5 KB

Line
1	!
2	! $Id: etat0_netcdf.F 1379 2010-05-06 12:19:18Z lguez $
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