source: LMDZ4/trunk/libf/dyn3dpar/etat0_netcdf.F90 @ 1319

Last change on this file since 1319 was 1319, checked in by Laurent Fairhead, 14 years ago
  • Modifications to the start and limit creation routines to account for different

calendars

  • Modification to phyetat0 to force the mask read in the start files to match the

surface fractions read in the limit file

  • Force readaerosol.F90 to read in aerosols file with 12 timesteps

  • Modifications aux routines de créations des fichiers start et limit pour prendre

en compte différents calendriers

  • Modification à phyetat0 pour forcer le masque lu dans le fichier start à être

compatible avec les fractions de surface lu dans le fichier limit

  • Forcer readaerosol à ne lire que des fichiers à 12 pas de temps
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 19.0 KB
Line 
1!
2! $Id: etat0_netcdf.F90 1319 2010-02-23 21:29:54Z fairhead $
3!
4!-------------------------------------------------------------------------------
5!
6SUBROUTINE 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(x, iip1, jjp1, rlonv, rlatu, masque, 0.0, jjm, rlonu,    &
174                  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    couple=.true.
182    iret=NF90_CLOSE(nid_o2a)
183    CALL flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp, nid_o2a)
184    IF(iml_omask/=iim .OR.jml_omask/=jjp1) THEN
185      WRITE(lunout,*)'Dimensions non compatibles pour masque ocean'
186      WRITE(lunout,*)'iim = ',iim,' iml_omask = ',iml_omask
187      WRITE(lunout,*)'jjp1 = ',jjp1,' jml_omask = ',jml_omask
188      CALL abort_gcm('etat0_netcdf','Dimensions non compatibles pour masque oc&
189     &ean',1)
190    END IF
191    ALLOCATE(  ocemask(iml_omask,jml_omask),   ocetmp(iml_omask,jml_omask))
192    ALLOCATE(lon_omask(iml_omask,jml_omask),lat_omask(iml_omask,jml_omask))
193    ALLOCATE(dlon_omask(iml_omask),dlat_omask(jml_omask))
194    CALL flinopen("o2a.nc", .FALSE., iml_omask, jml_omask, llm_tmp,            &
195                   lon_omask, lat_omask, lev, ttm_tmp, itaul, date, dt, fid)
196    CALL flinget(fid,'OceMask',iml_omask,jml_omask,llm_tmp,ttm_tmp,1,1,ocetmp)
197    CALL flinclo(fid)
198    dlon_omask(:)=lon_omask(:,1)
199    dlat_omask(:)=lat_omask(1,:)
200    ocemask=ocetmp
201    IF(dlat_omask(1)<dlat_omask(jml_omask)) THEN
202      DO j=1,jml_omask
203        ocemask(:,j) = ocetmp(:,jml_omask-j+1)
204      END DO
205    END IF
206    ALLOCATE(   ocemask(iml_omask,jml_omask),   ocetmp(iml_omask,jml_omask))
207    ALLOCATE( lon_omask(iml_omask,jml_omask),lat_omask(iml_omask,jml_omask))
208    ALLOCATE(dlon_omask(iml_omask),         dlat_omask(jml_omask))
209    CALL flinopen("o2a.nc", .FALSE., iml_omask, jml_omask, llm_tmp, lon_omask, &
210                  lat_omask, lev, ttm_tmp, itaul, date, dt, fid)
211    CALL flinget(fid, 'OceMask', iml_omask, jml_omask, llm_tmp, ttm_tmp,       &
212                  1, 1, ocetmp)
213    CALL flinclo(fid)
214    dlon_omask(1:iml_omask) = lon_omask(1:iml_omask,1)
215    dlat_omask(1:jml_omask) = lat_omask(1,1:jml_omask)
216    ocemask = ocetmp
217    IF(dlat_omask(1)<dlat_omask(jml_omask)) THEN
218      DO j=1,jml_omask
219        ocemask(:,j) = ocetmp(:,jml_omask-j+1)
220      END DO
221    END IF
222!
223! Ocean mask to physical grid
224!*******************************************************************************
225    WRITE(lunout,*)'ocemask '
226    WRITE(fmt,"(i4,'i1)')")iml_omask ; fmt='('//ADJUSTL(fmt)
227    WRITE(lunout,fmt)int(ocemask)
228    ocemask_fi(1)=ocemask(1,1)
229    DO j=2,jjm; ocemask_fi((j-2)*iim+1:iim+1)=ocemask(1:iim,j); END DO
230    ocemask_fi(klon)=ocemask(1,jjp1)
231    zmasq=1.-ocemask_fi
232  END IF
233
234  CALL gr_fi_dyn(1,klon,iip1,jjp1,zmasq,masque)
235
236  ! The startget calls need to be replaced by a call to restget to get the
237  ! values in the restart file
238  x = 'relief'; orog(:,:) = 0.0
239  CALL startget(x,iip1,jjp1,rlonv,rlatu,          orog, 0.0,jjm,rlonu,rlatv,ib,&
240                masque)
241
242  x = 'rugosite'; rugo(:,:) = 0.0
243  CALL startget(x,iip1,jjp1,rlonv,rlatu,          rugo, 0.0,jjm, rlonu,rlatv,ib)
244!  WRITE(lunout,'(49I1)') INT(orog(:,:)*10)
245!  WRITE(lunout,'(49I1)') INT(rugo(:,:)*10)
246
247! Sub-surfaces initialization
248!*******************************************************************************
249  pctsrf=0.
250  x = 'psol'; psol(:,:) = 0.0
251  CALL startget(x,iip1,jjp1,rlonv,rlatu,psol,0.0,jjm,rlonu,rlatv,ib)
252!  WRITE(lunout,*) 'PSOL :', psol(10,20)
253!  WRITE(lunout,*) ap(:), bp(:)
254
255! Mid-levels pressure computation
256!*******************************************************************************
257  CALL pression(ip1jmp1, ap, bp, psol, p3d)
258  CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, y)
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(x,iip1,jjp1,rlonv,rlatu,          phis, 0.0,jjm, rlonu,rlatv,ib)
266
267  x = 'u';    uvent(:,:,:) = 0.0
268  CALL startget(x,iip1,jjp1,rlonu,rlatu,llm,pls,y,uvent,0.0,jjm, rlonv,rlatv,ib)
269
270  x = 'v';    vvent(:,:,:) = 0.0
271  CALL startget(x,iip1,jjm, rlonv,rlatv,llm,pls,y,vvent,0.0,jjp1,rlonu,rlatu,ib)
272
273  x = 't';    t3d(:,:,:) = 0.0
274  CALL startget(x,iip1,jjp1,rlonv,rlatu,llm,pls,y,t3d,  0.0,jjm, rlonu,rlatv,ib)
275
276  x = 'tpot'; tpot(:,:,:) = 0.0
277  CALL startget(x,iip1,jjp1,rlonv,rlatu,llm,pls,pk,tpot,0.0,jjm, rlonu,rlatv,ib)
278
279  WRITE(lunout,*) 'T3D min,max:',minval(t3d(:,:,:)),maxval(t3d(:,:,:))
280  WRITE(lunout,*) 'PLS min,max:',minval(pls(:,:,:)),maxval(pls(:,:,:))
281
282! Humidity at saturation computation
283!*******************************************************************************
284  WRITE(lunout,*) 'avant q_sat'
285  CALL q_sat(llm*jjp1*iip1, t3d, pls, qsat)
286  WRITE(lunout,*) 'apres q_sat'
287  WRITE(lunout,*) 'QSAT min,max:',minval(qsat(:,:,:)),maxval(qsat(:,:,:))
288!  WRITE(lunout,*) 'QSAT :',qsat(10,20,:)
289
290  x = 'q';    qd (:,:,:) = 0.0
291  CALL startget(x,iip1,jjp1,rlonv,rlatu,llm,pls,qsat,qd,0.0,jjm, rlonu,rlatv,ib)
292  q3d(:,:,:,:) = 0.0 ; q3d(:,:,:,1) = qd(:,:,:)
293
294!--- OZONE CLIMATOLOGY
295  IF(read_climoz>=1) CALL regr_lat_time_climoz(read_climoz)
296
297  x = 'tsol'; tsol(:) = 0.0
298  CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,tsol,  0.0,jjm,rlonu,rlatv,ib)
299
300  x = 'qsol';  qsol(:) = 0.0
301  CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,qsol,  0.0,jjm,rlonu,rlatv,ib)
302
303  x = 'snow';  sn(:) = 0.0
304  CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,sn,    0.0,jjm,rlonu,rlatv,ib)
305
306  x = 'rads';  radsol(:) = 0.0
307  CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,radsol,0.0,jjm,rlonu,rlatv,ib)
308
309  x = 'rugmer'; rugmer(:) = 0.0
310  CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,rugmer,0.0,jjm,rlonu,rlatv,ib)
311
312  x = 'zmea';  zmea(:) = 0.0
313  CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,zmea,  0.0,jjm,rlonu,rlatv,ib)
314
315  x = 'zstd';  zstd(:) = 0.0
316  CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,zstd,  0.0,jjm,rlonu,rlatv,ib)
317
318  x = 'zsig';  zsig(:) = 0.0
319  CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,zsig,  0.0,jjm,rlonu,rlatv,ib)
320
321  x = 'zgam';  zgam(:) = 0.0
322  CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,zgam,  0.0,jjm,rlonu,rlatv,ib)
323
324  x = 'zthe';  zthe(:) = 0.0
325  CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,zthe,  0.0,jjm,rlonu,rlatv,ib)
326
327  x = 'zpic';  zpic(:) = 0.0
328  CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,zpic,  0.0,jjm,rlonu,rlatv,ib)
329
330  x = 'zval';  zval(:) = 0.0
331  CALL startget(x,iip1,jjp1,rlonv,rlatu,klon,zval,  0.0,jjm,rlonu,rlatv,ib)
332
333!  WRITE(lunout,'(48I3)') 'TSOL :', INT(tsol(2:klon)-273)
334
335! Soil ice file reading for soil fraction and soil ice fraction
336!*******************************************************************************
337  CALL flininfo("landiceref.nc", iml_lic, jml_lic, llm_tmp, ttm_tmp, fid)
338  ALLOCATE( lat_lic(iml_lic,jml_lic),lon_lic(iml_lic, jml_lic))
339  ALLOCATE(dlat_lic(jml_lic),       dlon_lic(iml_lic))
340  ALLOCATE( fraclic(iml_lic,jml_lic))
341  CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp,           &
342                lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid)
343  CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp, 1,1, fraclic)
344  CALL flinclo(fid)
345
346! Interpolation on model T-grid
347!*******************************************************************************
348  WRITE(lunout,*)'dimensions de landice iml_lic, jml_lic : ',iml_lic,jml_lic
349! conversion if coordinates are in degrees
350  IF(MAXVAL(lon_lic)>pi) lon_lic=lon_lic*pi/180.
351  IF(MAXVAL(lat_lic)>pi) lat_lic=lat_lic*pi/180.
352  dlon_lic(:)=lon_lic(:,1)
353  dlat_lic(:)=lat_lic(1,:)
354  CALL grille_m( iml_lic, jml_lic, dlon_lic, dlat_lic, fraclic, iim,jjp1,      &
355                 rlonv, rlatu, flic_tmp(1:iim,:) )
356  flic_tmp(iip1,:)=flic_tmp(1,:)
357
358!--- To the physical grid
359  CALL gr_dyn_fi(1, iip1, jjp1, klon, flic_tmp, pctsrf(:,is_lic))
360
361!--- Adequation with soil/sea mask
362  WHERE(pctsrf(:,is_lic)<EPSFRA) pctsrf(:,is_lic)=0.
363  WHERE(zmasq(:)<EPSFRA)         pctsrf(:,is_lic)=0.
364  pctsrf(:,is_ter)=zmasq(:)
365  DO ji=1,klon
366    IF(zmasq(ji)>EPSFRA) THEN
367      IF(pctsrf(ji,is_lic)>=zmasq(ji)) THEN
368        pctsrf(ji,is_lic)=zmasq(ji)
369        pctsrf(ji,is_ter)=0.
370      ELSE
371        pctsrf(ji,is_ter)=zmasq(ji)-pctsrf(ji,is_lic)
372        IF(pctsrf(ji,is_ter)<EPSFRA) THEN
373          pctsrf(ji,is_ter)=0.
374          pctsrf(ji,is_lic)=zmasq(ji)
375        END IF
376      END IF
377    END IF
378  END DO
379
380! sub-surface ocean and sea ice (sea ice set to zero for start)
381!*******************************************************************************
382  pctsrf(:,is_oce)=(1.-zmasq(:))
383  WHERE(pctsrf(:,is_oce)<EPSFRA) pctsrf(:,is_oce)=0.
384  IF(couple) pctsrf(:,is_oce)=ocemask_fi(:)
385  isst=0
386  WHERE(pctsrf(2:klon-1,is_oce)>0.) isst=1
387
388! It is checked that the sub-surfaces sum is equal to 1
389!*******************************************************************************
390  ji=COUNT((ABS(SUM(pctsrf(:,:),dim=2))-1.0)>EPSFRA)
391  IF(ji/=0) WRITE(lunout,*) 'pb repartition sous maille pour ',ji,' points'
392  CALL gr_fi_ecrit(1, klon, iim, jjp1, zmasq, zx_tmp_2d)
393!  WRITE(fmt,"(i3,')')")iim; fmt='(i'//ADJUSTL(fmt)
394!  WRITE(lunout,*)'zmasq = '
395!  WRITE(lunout,TRIM(fmt))NINT(zx_tmp_2d)
396  CALL gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque)
397  WRITE(fmt,"(i4,'i1)')")iip1 ; fmt='('//ADJUSTL(fmt)
398  WRITE(lunout,*) 'MASQUE construit : Masque'
399  WRITE(lunout,TRIM(fmt))NINT(masque(:,:))
400
401! Intermediate computation
402!*******************************************************************************
403  CALL massdair(p3d,masse)
404  WRITE(lunout,*)' ALPHAX ',alphax
405  DO l=1,llm
406    xppn(:)=aire(1:iim,1   )*masse(1:iim,1   ,l)
407    xpps(:)=aire(1:iim,jjp1)*masse(1:iim,jjp1,l)
408    xpn=SUM(xppn)/apoln
409    xps=SUM(xpps)/apols
410    masse(:,1   ,l)=xpn
411    masse(:,jjp1,l)=xps
412  END DO
413  q3d(iip1,:,:,:)=q3d(1,:,:,:)
414  phis(iip1,:) = phis(1,:)
415
416  IF(.NOT.letat0) RETURN
417
418! Writing
419!*******************************************************************************
420  CALL inidissip(lstardis,nitergdiv,nitergrot,niterh,tetagdiv,tetagrot,tetatemp)
421  WRITE(lunout,*)'sortie inidissip'
422  itau=0
423  itau_dyn=0
424  itau_phy=0
425  iday=dayref+itau/day_step
426  time=FLOAT(itau-(iday-dayref)*day_step)/day_step
427  IF(time>1.) THEN
428   time=time-1
429   iday=iday+1
430  END IF
431  day_ref=dayref
432  annee_ref=anneeref
433
434  CALL geopot( ip1jmp1, tpot, pk, pks, phis, phi )
435  WRITE(lunout,*)'sortie geopot'
436
437  CALL caldyn0( itau, uvent, vvent, tpot, psol, masse, pk, phis,               &
438                phi,  w, pbaru, pbarv, time+iday-dayref)
439  WRITE(lunout,*)'sortie caldyn0'     
440  CALL dynredem0( "start.nc", dayref, phis)
441  WRITE(lunout,*)'sortie dynredem0'
442  CALL dynredem1( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol)
443  WRITE(lunout,*)'sortie dynredem1'
444
445! Physical initial state writting
446!*******************************************************************************
447  WRITE(lunout,*)'phystep ',dtvr,iphysiq,nbapp_rad
448  phystep   = dtvr * FLOAT(iphysiq)
449  radpas    = NINT (86400./phystep/ FLOAT(nbapp_rad) )
450  WRITE(lunout,*)'phystep =', phystep, radpas
451
452! Init: tsol, qsol, sn, evap, tsoil, rain_fall, snow_fall, solsw, sollw, frugs
453!*******************************************************************************
454  DO i=1,nbsrf; ftsol(:,i) = tsol; END DO
455  DO i=1,nbsrf; snsrf(:,i) = sn;   END DO
456  falb1(:,is_ter) = 0.08; falb1(:,is_lic) = 0.6
457  falb1(:,is_oce) = 0.5;  falb1(:,is_sic) = 0.6
458  falb2 = falb1
459  evap(:,:) = 0.
460  DO i=1,nbsrf; qsolsrf(:,i)=150.; END DO
461  DO i=1,nbsrf; DO j=1,nsoilmx; tsoil(:,j,i) = tsol; END DO; END DO
462  rain_fall = 0.; snow_fall = 0.
463  solsw = 165.;   sollw = -53.
464  t_ancien = 273.15
465  q_ancien = 0.
466  agesno = 0.
467  frugs(:,is_oce) = rugmer(:)
468  frugs(:,is_ter) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
469  frugs(:,is_lic) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
470  frugs(:,is_sic) = 0.001
471  fder = 0.0
472  clwcon = 0.0
473  rnebcon = 0.0
474  ratqs = 0.0
475  run_off_lic_0 = 0.0
476  rugoro = 0.0
477
478! Before phyredem calling, surface modules and values to be saved in startphy.nc
479! are initialized
480!*******************************************************************************
481  dummy = 1.0
482  pbl_tke(:,:,:) = 1.e-8
483  zmax0(:) = 40.
484  f0(:) = 1.e-5
485  ema_work1(:,:) = 0.
486  ema_work2(:,:) = 0.
487  wake_deltat(:,:) = 0.
488  wake_deltaq(:,:) = 0.
489  wake_s(:) = 0.
490  wake_cstar(:) = 0.
491  wake_fip(:) = 0.
492  CALL fonte_neige_init(run_off_lic_0)
493  CALL pbl_surface_init( qsol, fder, snsrf, qsolsrf, evap, frugs, agesno, tsoil )
494  CALL phyredem( "startphy.nc" )
495
496  WRITE(lunout,*)'CCCCCCCCCCCCCCCCCC REACTIVER SORTIE VISU DANS ETAT0'
497  WRITE(lunout,*)'entree histclo'
498  CALL histclo()
499
500#endif
501!#endif of #ifdef CPP_EARTH
502  RETURN
503
504END SUBROUTINE etat0_netcdf
505!
506!-------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.