Context Navigation

etat0_netcdf.F90 @ 5072

Last change on this file since 5072 was 1520, checked in by Ehouarn Millour, 13 years ago

Implementation of a different vertical discretization (from/for planets, but
can in principle also be used for Earth).
Choice of vertical discretization is set by flag 'disvert_type';
'disvert_type=1' is Earth standard (default; ie set to 1 if
planet_type=="earth") case.
With 'disvert_type=2', approximate altitude of layers and reference atmospheric
scale height must be given using an input file ("z2sig.def", first line
should give scale height, in km, following lines must specify the altitude,
in km above surface, of mid-layers, one per line; see disvert_noterre.F).

Checked that these changes do not impact on 'bench' results, on Vargas.

EM.

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

File size: 19.4 KB

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

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

Download in other formats:

Original Format