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

Last change on this file since 2239 was 2239, checked in by Ehouarn Millour, 10 years ago

Reorganizing physics/dynamics interface:

  • what is related to dynamics-physics interface is now in a seperate directory: dynlmdz_phy* for physics in phy*
  • 1d model and related dependencies (including a couple from "dynamics", set up as symbolic links) is now in subdirectory "dyn1d" of phy*.
  • "bibio" directory is now "misc" and should only contain autonomous utilities.
  • "cosp" is now a subdirectory of phylmd.

EM

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