source: LMDZ4/branches/LMDZ4_AR5/libf/dyn3d/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!
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_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(llm*jjp1*iip1, 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
495END SUBROUTINE etat0_netcdf
496!
497!-------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.