source: LMDZ5/branches/LMDZ5_AR5/libf/dyn3d/etat0_netcdf.F90 @ 4116

Last change on this file since 4116 was 1520, checked in by Ehouarn Millour, 13 years ago

Implementation of a different vertical discretization (from/for planets, but
can in principle also be used for Earth).
Choice of vertical discretization is set by flag 'disvert_type';
'disvert_type=1' is Earth standard (default; ie set to 1 if
planet_type=="earth") case.
With 'disvert_type=2', approximate altitude of layers and reference atmospheric
scale height must be given using an input file ("z2sig.def", first line
should give scale height, in km, following lines must specify the altitude,
in km above surface, of mid-layers, one per line; see disvert_noterre.F).

Checked that these changes do not impact on 'bench' results, on Vargas.

EM.

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