source: LMDZ4/trunk/libf/dyn3d/etat0_netcdf.F90 @ 2916

Last change on this file since 2916 was 1425, checked in by lguez, 14 years ago

Replaced Numerical Recipes procedures for spline interpolation (not in
the public domain) by procedures from the Pchip package in the Slatec
library. This only affects the program "ce0l", not the program
"gcm". Tested on Brodie SX8 with "-debug" and "-prod", "-parallel
none" and "-parallel mpi". "start.nc" and "limit.nc" are
changed. "startphy.nc" is not changed. The relative change is of order
1e-7 or less. The revision makes the program faster (tested on Brodie
with "-prod -d 144x142x39", CPU time is 38 s, instead of 54
s). Procedures from Slatec are untouched, except for
"i1mach.F". Created procedures "pchfe_95" and "pchsp_95" which are
wrappers for "pchfe" and "pchsp" from Slatec. "pchfe_95" and
"pchsp_95" have a safer and simpler interface.

Replaced "make" by "sxgmake" in "arch-SX8_BRODIE.fcm". Added files for
compilation by FCM with "g95".

In "arch-linux-32bit.fcm", replaced "pgf90" by "pgf95". There was no
difference between "dev" and "debug" so added "-O1" to "dev". Added
debugging options. Removed "-Wl,-Bstatic
-L/usr/lib/gcc-lib/i386-linux/2.95.2", which usually produces an error
at link-time.

Bash is now ubiquitous while KornShell? is not so use Bash instead of
KornShell? in FCM.

Replaced some statements "write(6,*)" by "write(lunout,*)". Replaced
"stop" by "stop 1" in the case where "abort_gcm" is called with "ierr
/= 0". Removed "stop" statements at the end of procedures
"limit_netcdf" and main program "ce0l" (why not let the program end
normally?).

Made some arrays automatic instead of allocatable in "start_inter_3d".

Zeroed "wake_pe", "fm_therm", "entr_therm" and "detr_therm" in
"dyn3dpar/etat0_netcdf.F90". The parallel and sequential results of
"ce0l" are thus identical.

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