source: LMDZ5/branches/testing/libf/phylmd/etat0_netcdf.F90 @ 2056

Last change on this file since 2056 was 2056, checked in by Laurent Fairhead, 10 years ago

Merged trunk changes r1997:2055 into testing branch

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