Context Navigation

etat0_netcdf.F90 @ 5017

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

Continent/ocean mask was transformed incorrectly from o2a.nc file
Change of unit for seaice depending on whether the initial data is from
climatology or model results

Erreur sur la transformation vers la grille physique du masque lu dans o2a.nc
Pas de division par 100 de la glace océanique si on lit des champs du modèle
couplé

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

File size: 18.5 KB

Line
1	!
2	! $Id: etat0_netcdf.F90 1328 2010-03-18 13:26:23Z dcugnet $
3	!
4	!-------------------------------------------------------------------------------
5	!
6	SUBROUTINE etat0_netcdf(ib, masque, letat0)
7	!
8	!-------------------------------------------------------------------------------
9	! Purpose: Creates initial states
10	!-------------------------------------------------------------------------------
11	! Note: This routine is designed to work for Earth
12	!-------------------------------------------------------------------------------
13	#ifdef CPP_EARTH
14	USE startvar
15	USE ioipsl
16	USE dimphy
17	USE infotrac
18	USE fonte_neige_mod
19	USE pbl_surface_mod
20	USE phys_state_var_mod
21	USE filtreg_mod
22	USE regr_lat_time_climoz_m, ONLY: regr_lat_time_climoz
23	USE conf_phys_m, ONLY: conf_phys
24	USE netcdf, ONLY : NF90_OPEN, NF90_NOWRITE, NF90_CLOSE, NF90_NOERR
25	#endif
26	IMPLICIT NONE
27	!-------------------------------------------------------------------------------
28	! Arguments:
29	#include "dimensions.h"
30	#include "paramet.h"
31	#include "iniprint.h"
32	LOGICAL, INTENT(IN) :: ib ! barycentric interpolat.
33	REAL, DIMENSION(iip1,jjp1), INTENT(INOUT) :: masque ! land mask
34	LOGICAL, INTENT(IN) :: letat0 ! F: masque only required
35	#ifndef CPP_EARTH
36	WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics'
37	#else
38	!-------------------------------------------------------------------------------
39	! Local variables:
40	#include "comgeom2.h"
41	#include "comvert.h"
42	#include "comconst.h"
43	#include "indicesol.h"
44	#include "dimsoil.h"
45	#include "temps.h"
46	REAL, DIMENSION(klon) :: tsol, qsol
47	REAL, DIMENSION(klon) :: sn, rugmer, run_off_lic_0
48	REAL, DIMENSION(iip1,jjp1) :: orog, rugo, psol, phis
49	REAL, DIMENSION(iip1,jjp1,llm+1) :: p3d
50	REAL, DIMENSION(iip1,jjp1,llm) :: uvent, t3d, tpot, qsat, qd
51	REAL, DIMENSION(iip1,jjm ,llm) :: vvent
52	REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: q3d
53	REAL, DIMENSION(klon,nbsrf) :: qsolsrf, snsrf, evap
54	REAL, DIMENSION(klon,nbsrf) :: frugs, agesno
55	REAL, DIMENSION(klon,nsoilmx,nbsrf) :: tsoil
56
57	!--- Local variables for sea-ice reading:
58	INTEGER :: iml_lic, jml_lic, llm_tmp
59	INTEGER :: ttm_tmp, iret, fid
60	INTEGER, DIMENSION(1) :: itaul
61	REAL, DIMENSION(1) :: lev
62	REAL :: date
63	REAL, DIMENSION(:,:), ALLOCATABLE :: lon_lic, lat_lic, fraclic
64	REAL, DIMENSION(:), ALLOCATABLE :: dlon_lic, dlat_lic
65	REAL, DIMENSION(iip1,jjp1) :: flic_tmp
66
67	!--- Misc
68	CHARACTER(LEN=80) :: x, fmt
69	INTEGER :: i, j, l, ji
70	REAL, DIMENSION(iip1,jjp1,llm) :: alpha, beta, pk, pls, y
71	REAL, DIMENSION(ip1jmp1) :: pks
72
73	#include "comdissnew.h"
74	#include "control.h"
75	#include "serre.h"
76	#include "clesphys.h"
77
78	REAL, DIMENSION(iip1,jjp1,llm) :: masse
79	INTEGER :: itau, iday
80	REAL :: xpn, xps, time, phystep
81	REAL, DIMENSION(iim) :: xppn, xpps
82	REAL, DIMENSION(ip1jmp1,llm) :: pbaru, phi, w
83	REAL, DIMENSION(ip1jm ,llm) :: pbarv
84	REAL, DIMENSION(klon) :: fder
85
86	!--- Local variables for ocean mask reading:
87	INTEGER :: nid_o2a, iml_omask, jml_omask
88	LOGICAL :: couple=.FALSE.
89	REAL, DIMENSION(:,:), ALLOCATABLE :: lon_omask, lat_omask, ocemask, ocetmp
90	REAL, DIMENSION(:), ALLOCATABLE :: dlon_omask,dlat_omask
91	REAL, DIMENSION(klon) :: ocemask_fi
92	INTEGER, DIMENSION(klon-2) :: isst
93	REAL, DIMENSION(iim,jjp1) :: zx_tmp_2d
94	REAL :: dummy
95	LOGICAL :: ok_newmicro, ok_journe, ok_mensuel, ok_instan, ok_hf
96	LOGICAL :: ok_LES, ok_ade, ok_aie, aerosol_couple, new_aod
97	INTEGER :: iflag_radia, flag_aerosol
98	REAL :: bl95_b0, bl95_b1, fact_cldcon, facttemps, ratqsbas, ratqshaut
99	REAL :: tau_ratqs
100	INTEGER :: iflag_cldcon, iflag_ratqs, iflag_coupl, iflag_clos, iflag_wake
101	INTEGER :: iflag_thermals, nsplit_thermals
102	INTEGER :: iflag_thermals_ed, iflag_thermals_optflux
103	REAL :: tau_thermals, solarlong0, seuil_inversion
104	INTEGER :: read_climoz ! read ozone climatology
105	! Allowed values are 0, 1 and 2
106	! 0: do not read an ozone climatology
107	! 1: read a single ozone climatology that will be used day and night
108	! 2: read two ozone climatologies, the average day and night
109	! climatology and the daylight climatology
110	!-------------------------------------------------------------------------------
111	!--- Constants
112	pi = 4. * ATAN(1.)
113	rad = 6371229.
114	daysec = 86400.
115	omeg = 2.*pi/daysec
116	g = 9.8
117	kappa = 0.2857143
118	cpp = 1004.70885
119	preff = 101325.
120	pa = 50000.
121	jmp1 = jjm + 1
122
123	!--- CONSTRUCT A GRID
124	CALL conf_phys( ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES, &
125	solarlong0,seuil_inversion, &
126	fact_cldcon, facttemps,ok_newmicro,iflag_radia, &
127	iflag_cldcon, &
128	iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
129	ok_ade, ok_aie, aerosol_couple, &
130	flag_aerosol, new_aod, &
131	bl95_b0, bl95_b1, &
132	iflag_thermals,nsplit_thermals,tau_thermals, &
133	iflag_thermals_ed,iflag_thermals_optflux, &
134	iflag_coupl,iflag_clos,iflag_wake, read_climoz )
135
136	! co2_ppm0 : initial value of atmospheric CO2 from .def file (co2_ppm value)
137	co2_ppm0 = co2_ppm
138
139	dtvr = daysec/FLOAT(day_step)
140	WRITE(lunout,*)'dtvr',dtvr
141
142	CALL iniconst()
143	CALL inigeom()
144
145	!--- Initializations for tracers
146	CALL infotrac_init
147	ALLOCATE(q3d(iip1,jjp1,llm,nqtot))
148
149	CALL inifilr()
150	CALL phys_state_var_init(read_climoz)
151
152	rlat(1) = ASIN(1.0)
153	DO j=2,jjm; rlat((j-2)iim+2:(j-1)iim+1)=rlatu(j); END DO
154	rlat(klon) = - ASIN(1.0)
155	rlat(:)=rlat(:)*(180.0/pi)
156
157	rlon(1) = 0.0
158	DO j=2,jjm; rlon((j-2)iim+2:(j-1)iim+1)=rlonv(1:iim); END DO
159	rlon(klon) = 0.0
160	rlon(:)=rlon(:)*(180.0/pi)
161
162	! For a coupled simulation, the ocean mask from ocean model is used to compute
163	! the weights an to insure ocean fractions are the same for atmosphere and ocean
164	! Otherwise, mask is created using Relief file.
165
166	WRITE(lunout,*)'Essai de lecture masque ocean'
167	iret = NF90_OPEN("o2a.nc", NF90_NOWRITE, nid_o2a)
168	IF(iret/=NF90_NOERR) THEN
169	WRITE(lunout,*)'ATTENTION!! pas de fichier o2a.nc trouve'
170	WRITE(lunout,*)'Run force'
171	x='masque'
172	masque(:,:)=0.0
173	CALL startget_phys2d(x, iip1, jjp1, rlonv, rlatu, masque, 0.0, jjm, &
174	& rlonu, rlatv, ib)
175	WRITE(lunout,*)'MASQUE construit : Masque'
176	WRITE(lunout,'(97I1)') nINT(masque)
177	CALL gr_dyn_fi(1, iip1, jjp1, klon, masque, zmasq)
178	WHERE( zmasq(:)<EPSFRA) zmasq(:)=0.
179	WHERE(1.-zmasq(:)<EPSFRA) zmasq(:)=1.
180	ELSE
181	WRITE(lunout,*)'ATTENTION!! fichier o2a.nc trouve'
182	WRITE(lunout,*)'Run couple'
183	couple=.true.
184	iret=NF90_CLOSE(nid_o2a)
185	CALL flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp, nid_o2a)
186	IF(iml_omask/=iim .OR.jml_omask/=jjp1) THEN
187	WRITE(lunout,*)'Dimensions non compatibles pour masque ocean'
188	WRITE(lunout,*)'iim = ',iim,' iml_omask = ',iml_omask
189	WRITE(lunout,*)'jjp1 = ',jjp1,' jml_omask = ',jml_omask
190	CALL abort_gcm('etat0_netcdf','Dimensions non compatibles pour masque oc&
191	&ean',1)
192	END IF
193	ALLOCATE( ocemask(iml_omask,jml_omask), ocetmp(iml_omask,jml_omask))
194	ALLOCATE( lon_omask(iml_omask,jml_omask),lat_omask(iml_omask,jml_omask))
195	ALLOCATE(dlon_omask(iml_omask), dlat_omask(jml_omask))
196	CALL flinopen("o2a.nc", .FALSE., iml_omask, jml_omask, llm_tmp, lon_omask,&
197	& lat_omask, lev, ttm_tmp, itaul, date, dt, fid)
198	CALL flinget(fid, 'OceMask', iml_omask, jml_omask, llm_tmp, ttm_tmp, &
199	& 1, 1, ocetmp)
200	CALL flinclo(fid)
201	dlon_omask(1:iml_omask) = lon_omask(1:iml_omask,1)
202	dlat_omask(1:jml_omask) = lat_omask(1,1:jml_omask)
203	ocemask = ocetmp
204	IF(dlat_omask(1)<dlat_omask(jml_omask)) THEN
205	DO j=1,jml_omask
206	ocemask(:,j) = ocetmp(:,jml_omask-j+1)
207	END DO
208	END IF
209	!
210	! Ocean mask to physical grid
211	!*******************************************************************************
212	WRITE(lunout,*)'ocemask '
213	WRITE(fmt,"(i4,'i1)')")iml_omask ; fmt='('//ADJUSTL(fmt)
214	WRITE(lunout,fmt)int(ocemask)
215	ocemask_fi(1)=ocemask(1,1)
216	DO j=2,jjm; ocemask_fi((j-2)iim+2:(j-1)iim+1)=ocemask(1:iim,j); END DO
217	ocemask_fi(klon)=ocemask(1,jjp1)
218	zmasq=1.-ocemask_fi
219	END IF
220
221	CALL gr_fi_dyn(1,klon,iip1,jjp1,zmasq,masque)
222
223	! The startget calls need to be replaced by a call to restget to get the
224	! values in the restart file
225	x = 'relief'; orog(:,:) = 0.0
226	CALL startget_phys2d(x,iip1,jjp1,rlonv,rlatu, orog, 0.0,jjm,rlonu,rlatv,ib,&
227	& masque)
228
229	x = 'rugosite'; rugo(:,:) = 0.0
230	CALL startget_phys2d(x,iip1,jjp1,rlonv,rlatu, rugo, 0.0,jjm, rlonu,rlatv,ib)
231	! WRITE(lunout,'(49I1)') INT(orog(:,:)*10)
232	! WRITE(lunout,'(49I1)') INT(rugo(:,:)*10)
233
234	! Sub-surfaces initialization
235	!*******************************************************************************
236	pctsrf=0.
237	x = 'psol'; psol(:,:) = 0.0
238	CALL startget_phys2d(x,iip1,jjp1,rlonv,rlatu,psol,0.0,jjm,rlonu,rlatv,ib)
239	! WRITE(lunout,*) 'PSOL :', psol(10,20)
240	! WRITE(lunout,*) ap(:), bp(:)
241
242	! Mid-levels pressure computation
243	!*******************************************************************************
244	CALL pression(ip1jmp1, ap, bp, psol, p3d)
245	CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, y)
246	pls(:,:,:)=preff(pk(:,:,:)/cpp)*(1./kappa)
247	! WRITE(lunout,*) 'P3D :', p3d(10,20,:)
248	! WRITE(lunout,*) 'PK:', pk(10,20,:)
249	! WRITE(lunout,*) 'PLS :', pls(10,20,:)
250
251	x = 'surfgeo'; phis(:,:) = 0.0
252	CALL startget_phys2d(x,iip1,jjp1,rlonv,rlatu,phis, 0.0,jjm, rlonu,rlatv,ib)
253
254	x = 'u'; uvent(:,:,:) = 0.0
255	CALL startget_dyn(x,rlonu,rlatu,pls,y,uvent,0.0, &
256	& rlonv,rlatv,ib)
257
258	x = 'v'; vvent(:,:,:) = 0.0
259	CALL startget_dyn(x, rlonv,rlatv,pls(:, :jjm, :),y(:, :jjm, :),vvent,0.0, &
260	& rlonu,rlatu(:jjm),ib)
261
262	x = 't'; t3d(:,:,:) = 0.0
263	CALL startget_dyn(x,rlonv,rlatu,pls,y,t3d,0.0, &
264	& rlonu,rlatv,ib)
265
266	x = 'tpot'; tpot(:,:,:) = 0.0
267	CALL startget_dyn(x,rlonv,rlatu,pls,pk,tpot,0.0, &
268	& rlonu,rlatv,ib)
269
270	WRITE(lunout,*) 'T3D min,max:',minval(t3d(:,:,:)),maxval(t3d(:,:,:))
271	WRITE(lunout,*) 'PLS min,max:',minval(pls(:,:,:)),maxval(pls(:,:,:))
272
273	! Humidity at saturation computation
274	!*******************************************************************************
275	WRITE(lunout,*) 'avant q_sat'
276	CALL q_sat(llmjjp1iip1, t3d, pls, qsat)
277	WRITE(lunout,*) 'apres q_sat'
278	WRITE(lunout,*) 'QSAT min,max:',minval(qsat(:,:,:)),maxval(qsat(:,:,:))
279	! WRITE(lunout,*) 'QSAT :',qsat(10,20,:)
280
281	x = 'q'; qd (:,:,:) = 0.0
282	CALL startget_dyn(x,rlonv,rlatu,pls,qsat,qd,0.0, rlonu,rlatv,ib)
283	q3d(:,:,:,:) = 0.0 ; q3d(:,:,:,1) = qd(:,:,:)
284
285	!--- OZONE CLIMATOLOGY
286	IF(read_climoz>=1) CALL regr_lat_time_climoz(read_climoz)
287
288	x = 'tsol'; tsol(:) = 0.0
289	CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,tsol,0.0,jjm,rlonu,rlatv,ib)
290
291	x = 'qsol'; qsol(:) = 0.0
292	CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,qsol,0.0,jjm,rlonu,rlatv,ib)
293
294	x = 'snow'; sn(:) = 0.0
295	CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,sn,0.0,jjm,rlonu,rlatv,ib)
296
297	x = 'rads'; radsol(:) = 0.0
298	CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,radsol,0.0,jjm,rlonu,rlatv,ib)
299
300	x = 'rugmer'; rugmer(:) = 0.0
301	CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,rugmer,0.0,jjm,rlonu,rlatv,ib)
302
303	x = 'zmea'; zmea(:) = 0.0
304	CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,zmea,0.0,jjm,rlonu,rlatv,ib)
305
306	x = 'zstd'; zstd(:) = 0.0
307	CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,zstd,0.0,jjm,rlonu,rlatv,ib)
308
309	x = 'zsig'; zsig(:) = 0.0
310	CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,zsig,0.0,jjm,rlonu,rlatv,ib)
311
312	x = 'zgam'; zgam(:) = 0.0
313	CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,zgam,0.0,jjm,rlonu,rlatv,ib)
314
315	x = 'zthe'; zthe(:) = 0.0
316	CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,zthe,0.0,jjm,rlonu,rlatv,ib)
317
318	x = 'zpic'; zpic(:) = 0.0
319	CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,zpic,0.0,jjm,rlonu,rlatv,ib)
320
321	x = 'zval'; zval(:) = 0.0
322	CALL startget_phys1d(x,iip1,jjp1,rlonv,rlatu,klon,zval,0.0,jjm,rlonu,rlatv,ib)
323
324	! WRITE(lunout,'(48I3)') 'TSOL :', INT(tsol(2:klon)-273)
325
326	! Soil ice file reading for soil fraction and soil ice fraction
327	!*******************************************************************************
328	CALL flininfo("landiceref.nc", iml_lic, jml_lic, llm_tmp, ttm_tmp, fid)
329	ALLOCATE( lat_lic(iml_lic,jml_lic),lon_lic(iml_lic, jml_lic))
330	ALLOCATE(dlat_lic(jml_lic), dlon_lic(iml_lic))
331	ALLOCATE( fraclic(iml_lic,jml_lic))
332	CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp, &
333	& lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid)
334	CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp, 1,1, fraclic)
335	CALL flinclo(fid)
336
337	! Interpolation on model T-grid
338	!*******************************************************************************
339	WRITE(lunout,*)'dimensions de landice iml_lic, jml_lic : ',iml_lic,jml_lic
340	! conversion if coordinates are in degrees
341	IF(MAXVAL(lon_lic)>pi) lon_lic=lon_lic*pi/180.
342	IF(MAXVAL(lat_lic)>pi) lat_lic=lat_lic*pi/180.
343	dlon_lic(:)=lon_lic(:,1)
344	dlat_lic(:)=lat_lic(1,:)
345	CALL grille_m( iml_lic, jml_lic, dlon_lic, dlat_lic, fraclic, iim,jjp1, &
346	& rlonv, rlatu, flic_tmp(1:iim,:) )
347	flic_tmp(iip1,:)=flic_tmp(1,:)
348
349	!--- To the physical grid
350	CALL gr_dyn_fi(1, iip1, jjp1, klon, flic_tmp, pctsrf(:,is_lic))
351
352	!--- Adequation with soil/sea mask
353	WHERE(pctsrf(:,is_lic)<EPSFRA) pctsrf(:,is_lic)=0.
354	WHERE(zmasq(:)<EPSFRA) pctsrf(:,is_lic)=0.
355	pctsrf(:,is_ter)=zmasq(:)
356	DO ji=1,klon
357	IF(zmasq(ji)>EPSFRA) THEN
358	IF(pctsrf(ji,is_lic)>=zmasq(ji)) THEN
359	pctsrf(ji,is_lic)=zmasq(ji)
360	pctsrf(ji,is_ter)=0.
361	ELSE
362	pctsrf(ji,is_ter)=zmasq(ji)-pctsrf(ji,is_lic)
363	IF(pctsrf(ji,is_ter)<EPSFRA) THEN
364	pctsrf(ji,is_ter)=0.
365	pctsrf(ji,is_lic)=zmasq(ji)
366	END IF
367	END IF
368	END IF
369	END DO
370
371	! sub-surface ocean and sea ice (sea ice set to zero for start)
372	!*******************************************************************************
373	pctsrf(:,is_oce)=(1.-zmasq(:))
374	WHERE(pctsrf(:,is_oce)<EPSFRA) pctsrf(:,is_oce)=0.
375	IF(couple) pctsrf(:,is_oce)=ocemask_fi(:)
376	isst=0
377	WHERE(pctsrf(2:klon-1,is_oce)>0.) isst=1
378
379	! It is checked that the sub-surfaces sum is equal to 1
380	!*******************************************************************************
381	ji=COUNT((ABS(SUM(pctsrf(:,:),dim=2))-1.0)>EPSFRA)
382	IF(ji/=0) WRITE(lunout,*) 'pb repartition sous maille pour ',ji,' points'
383	CALL gr_fi_ecrit(1, klon, iim, jjp1, zmasq, zx_tmp_2d)
384	! WRITE(fmt,"(i3,')')")iim; fmt='(i'//ADJUSTL(fmt)
385	! WRITE(lunout,*)'zmasq = '
386	! WRITE(lunout,TRIM(fmt))NINT(zx_tmp_2d)
387	CALL gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque)
388	WRITE(fmt,"(i4,'i1)')")iip1 ; fmt='('//ADJUSTL(fmt)
389	WRITE(lunout,*) 'MASQUE construit : Masque'
390	WRITE(lunout,TRIM(fmt))NINT(masque(:,:))
391
392	! Intermediate computation
393	!*******************************************************************************
394	CALL massdair(p3d,masse)
395	WRITE(lunout,*)' ALPHAX ',alphax
396	DO l=1,llm
397	xppn(:)=aire(1:iim,1 )*masse(1:iim,1 ,l)
398	xpps(:)=aire(1:iim,jjp1)*masse(1:iim,jjp1,l)
399	xpn=SUM(xppn)/apoln
400	xps=SUM(xpps)/apols
401	masse(:,1 ,l)=xpn
402	masse(:,jjp1,l)=xps
403	END DO
404	q3d(iip1,:,:,:)=q3d(1,:,:,:)
405	phis(iip1,:) = phis(1,:)
406
407	IF(.NOT.letat0) RETURN
408
409	! Writing
410	!*******************************************************************************
411	CALL inidissip(lstardis,nitergdiv,nitergrot,niterh,tetagdiv,tetagrot,tetatemp)
412	WRITE(lunout,*)'sortie inidissip'
413	itau=0
414	itau_dyn=0
415	itau_phy=0
416	iday=dayref+itau/day_step
417	time=FLOAT(itau-(iday-dayref)*day_step)/day_step
418	IF(time>1.) THEN
419	time=time-1
420	iday=iday+1
421	END IF
422	day_ref=dayref
423	annee_ref=anneeref
424
425	CALL geopot( ip1jmp1, tpot, pk, pks, phis, phi )
426	WRITE(lunout,*)'sortie geopot'
427
428	CALL caldyn0( itau, uvent, vvent, tpot, psol, masse, pk, phis, &
429	phi, w, pbaru, pbarv, time+iday-dayref)
430	WRITE(lunout,*)'sortie caldyn0'
431	CALL dynredem0( "start.nc", dayref, phis)
432	WRITE(lunout,*)'sortie dynredem0'
433	CALL dynredem1( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol)
434	WRITE(lunout,*)'sortie dynredem1'
435
436	! Physical initial state writting
437	!*******************************************************************************
438	WRITE(lunout,*)'phystep ',dtvr,iphysiq,nbapp_rad
439	phystep = dtvr * FLOAT(iphysiq)
440	radpas = NINT (86400./phystep/ FLOAT(nbapp_rad) )
441	WRITE(lunout,*)'phystep =', phystep, radpas
442
443	! Init: tsol, qsol, sn, evap, tsoil, rain_fall, snow_fall, solsw, sollw, frugs
444	!*******************************************************************************
445	DO i=1,nbsrf; ftsol(:,i) = tsol; END DO
446	DO i=1,nbsrf; snsrf(:,i) = sn; END DO
447	falb1(:,is_ter) = 0.08; falb1(:,is_lic) = 0.6
448	falb1(:,is_oce) = 0.5; falb1(:,is_sic) = 0.6
449	falb2 = falb1
450	evap(:,:) = 0.
451	DO i=1,nbsrf; qsolsrf(:,i)=150.; END DO
452	DO i=1,nbsrf; DO j=1,nsoilmx; tsoil(:,j,i) = tsol; END DO; END DO
453	rain_fall = 0.; snow_fall = 0.
454	solsw = 165.; sollw = -53.
455	t_ancien = 273.15
456	q_ancien = 0.
457	agesno = 0.
458	frugs(:,is_oce) = rugmer(:)
459	frugs(:,is_ter) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
460	frugs(:,is_lic) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
461	frugs(:,is_sic) = 0.001
462	fder = 0.0
463	clwcon = 0.0
464	rnebcon = 0.0
465	ratqs = 0.0
466	run_off_lic_0 = 0.0
467	rugoro = 0.0
468
469	! Before phyredem calling, surface modules and values to be saved in startphy.nc
470	! are initialized
471	!*******************************************************************************
472	dummy = 1.0
473	pbl_tke(:,:,:) = 1.e-8
474	zmax0(:) = 40.
475	f0(:) = 1.e-5
476	ema_work1(:,:) = 0.
477	ema_work2(:,:) = 0.
478	wake_deltat(:,:) = 0.
479	wake_deltaq(:,:) = 0.
480	wake_s(:) = 0.
481	wake_cstar(:) = 0.
482	wake_fip(:) = 0.
483	CALL fonte_neige_init(run_off_lic_0)
484	CALL pbl_surface_init( qsol, fder, snsrf, qsolsrf, evap, frugs, agesno, tsoil )
485	CALL phyredem( "startphy.nc" )
486
487	! WRITE(lunout,*)'CCCCCCCCCCCCCCCCCC REACTIVER SORTIE VISU DANS ETAT0'
488	! WRITE(lunout,*)'entree histclo'
489	CALL histclo()
490
491	#endif
492	!#endif of #ifdef CPP_EARTH
493	RETURN
494
495	END SUBROUTINE etat0_netcdf
496	!
497	!-------------------------------------------------------------------------------

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

Download in other formats:

Original Format