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

Last change on this file since 2337 was 2298, checked in by Laurent Fairhead, 9 years ago

Merged trunk changes -r2237:2291 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.6 KB
RevLine 
[1687]1!
2! $Id: etat0_netcdf.F90 1625 2012-05-09 13:14:48Z lguez $
3!
4!-------------------------------------------------------------------------------
5!
6SUBROUTINE etat0_netcdf(ib, masque, phis, letat0)
[2298]7#ifndef CPP_1D
[1687]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:
[1795]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
[1687]30  USE netcdf, ONLY : NF90_OPEN, NF90_NOWRITE, NF90_CLOSE, NF90_NOERR
[1795]31  USE indice_sol_mod
[2056]32  use exner_hyb_m, only: exner_hyb
33  use exner_milieu_m, only: exner_milieu
34  use test_disvert_m, only: test_disvert
[1687]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"
[2298]56  REAL,    DIMENSION(klon)                 :: tsol
[1687]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
[2298]63  REAL,    DIMENSION(klon,nbsrf)           :: qsolsrf, snsrf
[1687]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
[2056]79  REAL,    DIMENSION(iip1,jjp1,llm)        :: pk, pls, y
[1687]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
[1750]104  LOGICAL :: ok_LES, ok_ade, ok_aie, ok_cdnc, aerosol_couple, new_aod, callstats
[1687]105  INTEGER :: iflag_radia, flag_aerosol
[1795]106  LOGICAL :: flag_aerosol_strat
[1687]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,            &
[1750]142                   ok_ade, ok_aie, ok_cdnc, aerosol_couple,             &
[1795]143                   flag_aerosol, flag_aerosol_strat, new_aod,           &
[1687]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()
[2056]155  if (pressure_exner) call test_disvert
[1687]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'
[2160]189    if (prt_level >= 1) WRITE(lunout,'(97I1)') nINT(masque)
[1687]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
[2056]259    CALL exner_hyb(ip1jmp1, psol, p3d, pks, pk)
[1687]260  else
[2056]261    CALL exner_milieu(ip1jmp1,psol,p3d, pks,pk)
[1687]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'
[2160]423  if (prt_level >= 1) WRITE(lunout,TRIM(fmt))NINT(masque(:,:))
[1687]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!*******************************************************************************
[1750]444  CALL inidissip(lstardis, nitergdiv, nitergrot, niterh, tetagdiv, tetagrot, &
445       tetatemp, vert_prof_dissip)
[1687]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
[2258]481!albedo SB >>>
482  falb_dir(:,is_ter,:)=0.08; falb_dir(:,is_lic,:)=0.6
483  falb_dir(:,is_oce,:)=0.5;  falb_dir(:,is_sic,:)=0.6
484!albedo SB <<<
[2298]485  fevap(:,:) = 0.
[1687]486  DO i=1,nbsrf; qsolsrf(:,i)=150.; END DO
487  DO i=1,nbsrf; DO j=1,nsoilmx; tsoil(:,j,i) = tsol; END DO; END DO
488  rain_fall = 0.; snow_fall = 0.
489  solsw = 165.;   sollw = -53.
490  t_ancien = 273.15
491  q_ancien = 0.
492  agesno = 0.
[2298]493
494  z0m(:,is_oce) = rugmer(:)
495  z0m(:,is_ter) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
496  z0m(:,is_lic) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
497  z0m(:,is_sic) = 0.001
498  z0h(:,:)=z0m(:,:)
499
[1687]500  fder = 0.0
501  clwcon = 0.0
502  rnebcon = 0.0
503  ratqs = 0.0
504  run_off_lic_0 = 0.0
505  rugoro = 0.0
506
507! Before phyredem calling, surface modules and values to be saved in startphy.nc
508! are initialized
509!*******************************************************************************
510  dummy = 1.0
511  pbl_tke(:,:,:) = 1.e-8
512  zmax0(:) = 40.
513  f0(:) = 1.e-5
[1864]514  sig1(:,:) = 0.
515  w01(:,:) = 0.
[1687]516  wake_deltat(:,:) = 0.
517  wake_deltaq(:,:) = 0.
518  wake_s(:) = 0.
519  wake_cstar(:) = 0.
520  wake_fip(:) = 0.
521  wake_pe = 0.
522  fm_therm = 0.
523  entr_therm = 0.
524  detr_therm = 0.
525
526  CALL fonte_neige_init(run_off_lic_0)
[2298]527  CALL pbl_surface_init( fder, snsrf, qsolsrf, tsoil )
[1687]528  CALL phyredem( "startphy.nc" )
529
530!  WRITE(lunout,*)'CCCCCCCCCCCCCCCCCC REACTIVER SORTIE VISU DANS ETAT0'
531!  WRITE(lunout,*)'entree histclo'
532  CALL histclo()
533
534#endif
535!#endif of #ifdef CPP_EARTH
536  RETURN
[2298]537#endif
538!#endif of ifndef CPP_1D
[1687]539END SUBROUTINE etat0_netcdf
540!
541!-------------------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.