source: LMDZ5/trunk/libf/phylmd/etat0_netcdf.F90 @ 2019

Last change on this file since 2019 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

  • 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.4 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#endif
32  IMPLICIT NONE
33!-------------------------------------------------------------------------------
34! Arguments:
35#include "dimensions.h"
36#include "paramet.h"
37#include "iniprint.h"
38  LOGICAL,                    INTENT(IN)    :: ib     ! barycentric interpolat.
39  REAL, DIMENSION(iip1,jjp1), INTENT(INOUT) :: masque ! land mask
40  REAL, DIMENSION(iip1,jjp1), INTENT(OUT)   :: phis   ! geopotentiel au sol
41  LOGICAL,                    INTENT(IN)    :: letat0 ! F: masque only required
42#ifndef CPP_EARTH
43  WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics'
44#else
45!-------------------------------------------------------------------------------
46! Local variables:
47#include "comgeom2.h"
48#include "comvert.h"
49#include "comconst.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','Dimensions non compatibles pour masque oc&
200     &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  sig1(:,:) = 0.
507  w01(:,:) = 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
530END SUBROUTINE etat0_netcdf
531!
532!-------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.