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