source: LMDZ4/trunk/libf/dyn3d/etat0_netcdf.F @ 966

Last change on this file since 966 was 937, checked in by lmdzadmin, 17 years ago

Ajout variables convection (ema_work1, ema_work2) dans startphy.nc
IM
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 22.7 KB
Line 
1!
2! $Header$
3!
4c
5c
6      SUBROUTINE etat0_netcdf (interbar, masque, pctsrf)
7   
8      USE startvar
9      USE ioipsl
10      USE dimphy
11      USE fonte_neige_mod
12      USE pbl_surface_mod
13      USE ocean_slab_mod
14      USE surface_data
15      !
16      IMPLICIT NONE
17      !
18#include "netcdf.inc"
19#include "dimensions.h"
20#include "paramet.h"
21      !
22      !
23!      INTEGER, PARAMETER :: KIDIA=1, KFDIA=iim*(jjm-1)+2,
24!     .KLON=KFDIA-KIDIA+1,KLEV=llm
25      !
26#include "comgeom2.h"
27#include "comvert.h"
28#include "comconst.h"
29#include "indicesol.h"
30#include "dimsoil.h"
31#include "temps.h"
32      !
33      LOGICAL interbar
34      REAL :: latfi(klon), lonfi(klon)
35      REAL :: orog(iip1,jjp1), rugo(iip1,jjp1), masque(iip1,jjp1),
36     . psol(iip1, jjp1), phis(iip1, jjp1)
37      REAL :: p3d(iip1, jjp1, llm+1)
38      REAL :: uvent(iip1, jjp1, llm)
39      REAL :: vvent(iip1, jjm, llm)
40      REAL :: t3d(iip1, jjp1, llm), tpot(iip1, jjp1, llm)
41      REAL :: q3d(iip1, jjp1, llm,nqmx), qsat(iip1, jjp1, llm)
42      REAL :: tsol(klon), qsol(klon), sn(klon)
43      REAL :: tsolsrf(klon,nbsrf), qsolsrf(klon,nbsrf),snsrf(klon,nbsrf)
44      REAL :: albe(klon,nbsrf), evap(klon,nbsrf)
45      REAL :: alblw(klon,nbsrf)
46      REAL :: tsoil(klon,nsoilmx,nbsrf)
47      REAL :: radsol(klon),rain_fall(klon), snow_fall(klon)
48      REAL :: solsw(klon), sollw(klon), fder(klon)
49cIM "slab" ocean
50      REAL :: tslab(klon), seaice(klon)
51      REAL :: frugs(klon,nbsrf), agesno(klon,nbsrf)
52      REAL :: rugmer(klon)
53      REAL :: zmea(iip1*jjp1), zstd(iip1*jjp1)
54      REAL :: zsig(iip1*jjp1), zgam(iip1*jjp1), zthe(iip1*jjp1)
55      REAL :: zpic(iip1*jjp1), zval(iip1*jjp1), rugsrel(iip1*jjp1)
56      REAL :: qd(iip1, jjp1, llm)
57      REAL :: pctsrf(klon, nbsrf)
58      REAL :: t_ancien(klon,klev), q_ancien(klon,klev)      !
59      REAL :: run_off_lic_0(klon)
60      real :: clwcon(klon,klev),rnebcon(klon,klev),ratqs(klon,klev)
61      ! declarations pour lecture glace de mer
62      REAL :: rugv(klon)
63      INTEGER :: iml_lic, jml_lic, llm_tmp, ttm_tmp, iret
64      INTEGER :: itaul(1), fid
65      REAL :: lev(1), date
66      REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_lic, lat_lic
67      REAL, ALLOCATABLE, DIMENSION(:)  :: dlon_lic, dlat_lic
68      REAL, ALLOCATABLE, DIMENSION (:,:) :: fraclic
69      REAL :: flic_tmp(iip1, jjp1)
70      REAL :: champint(iim, jjp1)
71      !
72
73      CHARACTER*80 :: varname
74      !
75      INTEGER :: i,j, ig, l, ji,ii1,ii2
76      INTEGER :: nq
77      REAL :: xpi
78      !
79      REAL :: alpha(iip1,jjp1,llm),beta(iip1,jjp1,llm)
80      REAL :: pk(iip1,jjp1,llm), pls(iip1,jjp1,llm), pks(ip1jmp1)
81      REAL :: workvar(iip1,jjp1,llm)
82      !
83      REAL ::  prefkap, unskap
84      !
85      real :: time_step,t_ops,t_wrt
86
87#include "comdissnew.h"
88#include "control.h"
89#include "serre.h"
90#include "clesphys.h"
91
92      INTEGER  ::        longcles
93      PARAMETER      ( longcles  = 20 )
94      REAL :: clesphy0 ( longcles       )
95      REAL :: p(iip1,jjp1,llm)
96      INTEGER :: itau, iday
97      REAL :: masse(iip1,jjp1,llm)
98      REAL :: xpn,xps,xppn(iim),xpps(iim)
99      real :: time
100      REAL :: phi(ip1jmp1,llm)
101      REAL :: pbaru(ip1jmp1,llm),pbarv(ip1jm,llm)
102      REAL :: w(ip1jmp1,llm)
103      REAL ::phystep
104      INTEGER :: radpas
105       real zrel(iip1*jjp1),chmin,chmax
106
107      CHARACTER*80 :: visu_file
108      INTEGER :: visuid
109
110! pour la lecture du fichier masque ocean
111      integer :: nid_o2a
112      logical :: couple = .false.
113      INTEGER :: iml_omask, jml_omask
114      REAL, ALLOCATABLE, DIMENSION(:,:) :: lon_omask, lat_omask
115      REAL, ALLOCATABLE, DIMENSION(:)  :: dlon_omask, dlat_omask
116      REAL, ALLOCATABLE, DIMENSION (:,:) :: ocemask, ocetmp
117      real, dimension(klon) :: ocemask_fi
118      integer :: isst(klon-2)
119      real zx_tmp_2d(iim,jjp1)
120
121      REAL, DIMENSION(klon,klev,nbsrf) :: pbl_tke ! turb kinetic energy
122      REAL, DIMENSION(klon) :: zmax0, f0
123      REAL, DIMENSION(klon,klev) :: ema_work1, ema_work2
124      REAL :: dummy(klon, nbsrf)
125      !
126      !   Constantes
127      !
128      pi     = 4. * ATAN(1.)
129      rad    = 6371229.
130      omeg   = 4.* ASIN(1.)/(24.*3600.)
131      g      = 9.8
132      daysec = 86400.
133      kappa  = 0.2857143
134      cpp    = 1004.70885
135      !
136      preff     = 101325.
137      unskap = 1./kappa
138      !
139      jmp1    = jjm + 1
140      !
141      !    Construct a grid
142      !
143
144!      CALL defrun_new(99,.TRUE.,clesphy0)
145      CALL conf_gcm( 99, .TRUE. , clesphy0 )
146
147      dtvr   = daysec/FLOAT(day_step)
148      print*,'dtvr',dtvr
149
150      CALL inicons0()
151      CALL inigeom()
152      !
153      CALL inifilr()
154      !
155      latfi(1) = ASIN(1.0)
156      DO j = 2, jjm
157        DO i = 1, iim
158          latfi((j-2)*iim+1+i)=  rlatu(j)
159        ENDDO
160      ENDDO
161      latfi(klon) = - ASIN(1.0)
162      !
163      lonfi(1) = 0.0
164      DO j = 2, jjm
165        DO i = 1, iim
166          lonfi((j-2)*iim+1+i) =  rlonv(i)
167        ENDDO
168      ENDDO
169      lonfi(klon) = 0.0
170      !
171      xpi = 2.0 * ASIN(1.0)
172      DO ig = 1, klon
173        latfi(ig) = latfi(ig) * 180.0 / xpi
174        lonfi(ig) = lonfi(ig) * 180.0 / xpi
175      ENDDO
176      !
177
178
179C
180C En cas de simulation couplee, lecture du masque ocean issu du modele ocean
181C utilise pour calculer les poids et pour assurer l'adequation entre les
182C fractions d'ocean vu par l'atmosphere et l'ocean. Sinon, on cree le masque
183C a partir du fichier relief
184C
185
186      write(*,*)'Essai de lecture masque ocean'
187      iret = nf_open("o2a.nc", NF_NOWRITE, nid_o2a)
188      if (iret .ne. 0) then
189        write(*,*)'ATTENTION!! pas de fichier o2a.nc trouve'
190        write(*,*)'Run force'
191        varname = 'masque'
192        masque(:,:) = 0.0
193        CALL startget(varname, iip1, jjp1, rlonv, rlatu, masque, 0.0,
194     ,  jjm ,rlonu,rlatv , interbar )
195        WRITE(*,*) 'MASQUE construit : Masque'
196        WRITE(*,'(97I1)') nINT(masque(:,:))
197        call gr_dyn_fi(1, iip1, jjp1, klon, masque, zmasq)
198        WHERE (zmasq(1 : klon) .LT. EPSFRA)
199            zmasq(1 : klon) = 0.
200        END WHERE
201        WHERE (1. - zmasq(1 : klon) .LT. EPSFRA)
202            zmasq(1 : klon) = 1.
203        END WHERE
204      else
205        couple = .true.
206        iret = nf_close(nid_o2a)
207        call flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp
208     $    , nid_o2a)
209        if (iml_omask /= iim .or. jml_omask /= jjp1) then
210          write(*,*)'Dimensions non compatibles pour masque ocean'
211          write(*,*)'iim = ',iim,' iml_omask = ',iml_omask
212          write(*,*)'jjp1 = ',jjp1,' jml_omask = ',jml_omask
213          stop
214        endif
215        ALLOCATE(lat_omask(iml_omask, jml_omask), stat=iret)
216        ALLOCATE(lon_omask(iml_omask, jml_omask), stat=iret)
217        ALLOCATE(dlon_omask(iml_omask), stat=iret)
218        ALLOCATE(dlat_omask(jml_omask), stat=iret)
219        ALLOCATE(ocemask(iml_omask, jml_omask), stat=iret)
220        ALLOCATE(ocetmp(iml_omask, jml_omask), stat=iret)
221        CALL flinopen("o2a.nc", .FALSE., iml_omask, jml_omask, llm_tmp
222     $    , lon_omask, lat_omask, lev, ttm_tmp, itaul, date, dt, fid)
223        CALL flinget(fid, 'OceMask', iml_omask, jml_omask, llm_tmp,
224     $      ttm_tmp, 1, 1, ocetmp)
225        CALL flinclo(fid)
226        dlon_omask(1 : iml_omask) = lon_omask(1 : iml_omask, 1)
227        dlat_omask(1 : jml_omask) = lat_omask(1 , 1 : jml_omask)
228        ocemask = ocetmp
229        if (dlat_omask(1) < dlat_omask(jml_omask)) then
230          do j = 1, jml_omask
231            ocemask(:,j) = ocetmp(:,jml_omask-j+1)
232          enddo
233        endif
234C
235C passage masque ocean a la grille physique
236C
237        write(*,*)'ocemask '
238        write(*,'(96i1)')int(ocemask)
239        ocemask_fi(1) = ocemask(1,1)
240        do j = 2, jjm
241          do i = 1, iim
242            ocemask_fi((j-2)*iim + i + 1) = ocemask(i,j)
243          enddo
244        enddo
245        ocemask_fi(klon) = ocemask(1,jjp1)
246        zmasq = 1. - ocemask_fi
247      endif
248
249      call gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque)
250
251      varname = 'relief'
252      ! This line needs to be replaced by a call to restget to get the values in the restart file
253      orog(:,:) = 0.0
254       CALL startget(varname, iip1, jjp1, rlonv, rlatu, orog, 0.0 ,
255     , jjm ,rlonu,rlatv , interbar, masque )
256      !
257      WRITE(*,*) 'OUT OF GET VARIABLE : Relief'
258!      WRITE(*,'(49I1)') INT(orog(:,:))
259      !
260      varname = 'rugosite'
261      ! This line needs to be replaced by a call to restget to get the values in the restart file
262      rugo(:,:) = 0.0
263       CALL startget(varname, iip1, jjp1, rlonv, rlatu, rugo, 0.0 ,
264     , jjm, rlonu,rlatv , interbar )
265      !
266      WRITE(*,*) 'OUT OF GET VARIABLE : Rugosite'
267!      WRITE(*,'(49I1)') INT(rugo(:,:)*10)
268      !
269C
270C on initialise les sous surfaces
271C
272      pctsrf=0.
273c
274      varname = 'psol'
275      psol(:,:) = 0.0
276      CALL startget(varname, iip1, jjp1, rlonv, rlatu, psol, 0.0 ,
277     , jjm ,rlonu,rlatv , interbar )
278      !
279      !  Compute here the pressure on the intermediate levels. One would expect that this is available in the GCM
280      !  anyway.
281      !
282!      WRITE(*,*) 'PSOL :', psol(10,20)
283!      WRITE(*,*) ap(:), bp(:)
284      CALL pression(ip1jmp1, ap, bp, psol, p3d)
285!      WRITE(*,*) 'P3D :', p3d(10,20,:)
286      CALL exner_hyb(ip1jmp1, psol, p3d, alpha, beta, pks, pk, workvar)
287!      WRITE(*,*) 'PK:', pk(10,20,:)
288      !
289      !
290      !
291      prefkap =  preff  ** kappa
292!      WRITE(*,*) 'unskap, cpp,  preff :', unskap, cpp,  preff
293      DO l = 1, llm
294        DO j=1,jjp1
295          DO i =1, iip1
296            pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
297           ENDDO
298        ENDDO
299      ENDDO
300      !
301!      WRITE(*,*) 'PLS :', pls(10,20,:)
302      !
303      varname = 'surfgeo'
304      phis(:,:) = 0.0
305      CALL startget(varname, iip1, jjp1, rlonv, rlatu, phis, 0.0 ,
306     , jjm ,rlonu,rlatv, interbar )
307      !
308      varname = 'u'
309      uvent(:,:,:) = 0.0
310      CALL startget(varname, iip1, jjp1, rlonu, rlatu, llm, pls,
311     . workvar, uvent, 0.0, jjm ,rlonv, rlatv, interbar )
312      ! 
313      varname = 'v'
314      vvent(:,:,:) = 0.0
315      CALL startget(varname, iip1, jjm, rlonv, rlatv, llm, pls,
316     . workvar, vvent, 0.0, jjp1, rlonu, rlatu, interbar )
317      !
318      varname = 't'
319      t3d(:,:,:) = 0.0
320      CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
321     . workvar, t3d, 0.0 , jjm, rlonu, rlatv , interbar )
322      !
323      WRITE(*,*) 'T3D min,max:',minval(t3d(:,:,:)),
324     .                          maxval(t3d(:,:,:))
325      varname = 'tpot'
326      tpot(:,:,:) = 0.0
327      CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
328     . pk, tpot, 0.0 , jjm, rlonu, rlatv , interbar )
329      !
330      WRITE(*,*) 'T3D min,max:',minval(t3d(:,:,:)),
331     .                          maxval(t3d(:,:,:))
332      WRITE(*,*) 'PLS min,max:',minval(pls(:,:,:)),
333     .                          maxval(pls(:,:,:))
334
335c Calcul de l'humidite a saturation
336      print*,'avant q_sat'
337      call q_sat(llm*jjp1*iip1,t3d,pls,qsat)
338      print*,'apres q_sat'
339
340      WRITE(*,*) 'QSAT min,max:',minval(qsat(:,:,:)),
341     .                           maxval(qsat(:,:,:))
342      !
343      WRITE(*,*) 'QSAT :', qsat(10,20,:)
344      !
345      varname = 'q'
346      qd(:,:,:) = 0.0
347      q3d(:,:,:,:) = 0.0
348      WRITE(*,*) 'QSAT min,max:',minval(qsat(:,:,:)),
349     .                           maxval(qsat(:,:,:))
350      CALL startget(varname, iip1, jjp1, rlonv, rlatu, llm, pls,
351     . qsat, qd, 0.0, jjm, rlonu, rlatv , interbar )
352      q3d(:,:,:,1) = qd(:,:,:)
353      !
354      varname = 'tsol'
355      ! This line needs to be replaced by a call to restget to get the values in the restart file
356      tsol(:) = 0.0
357      CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, tsol, 0.0,
358     .    jjm, rlonu, rlatv , interbar )
359      !
360      WRITE(*,*) 'TSOL construit :'
361!      WRITE(*,'(48I3)') INT(TSOL(2:klon)-273)
362      !
363      varname = 'qsol'
364      qsol(:) = 0.0
365      CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, qsol, 0.0,
366     .   jjm, rlonu, rlatv , interbar )
367      !
368      varname = 'snow'
369      sn(:) = 0.0
370      CALL startget(varname, iip1, jjp1, rlonv, rlatu, klon, sn, 0.0,
371     .    jjm, rlonu, rlatv , interbar )
372      !
373      varname = 'rads'
374      radsol(:) = 0.0
375      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,radsol,0.0,
376     .    jjm, rlonu, rlatv , interbar )
377      !
378cIM "slab" ocean
379      varname = 'tslab'
380      tslab(:) = 0.0
381      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,tslab,0.0,
382     .     jjm, rlonu, rlatv , interbar )
383c
384      varname = 'seaice'
385      seaice(:) = 0.0
386      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,seaice,0.0,
387     .     jjm, rlonu, rlatv , interbar )
388      !
389      varname = 'rugmer'
390      rugmer(:) = 0.0
391      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,rugmer,0.0,
392     .     jjm, rlonu, rlatv , interbar )
393      !
394!      varname = 'agesno'
395!      agesno(:) = 0.0
396!      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,agesno,0.0,
397!     .     jjm, rlonu, rlatv , interbar )
398
399      varname = 'zmea'
400      zmea(:) = 0.0
401      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zmea,0.0,
402     .     jjm, rlonu, rlatv , interbar )
403
404      varname = 'zstd'
405      zstd(:) = 0.0
406      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zstd,0.0,
407     .     jjm, rlonu, rlatv , interbar )
408      varname = 'zsig'
409      zsig(:) = 0.0
410      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zsig,0.0,
411     .     jjm, rlonu, rlatv , interbar )
412      varname = 'zgam'
413      zgam(:) = 0.0
414      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zgam,0.0,
415     .     jjm, rlonu, rlatv , interbar )
416      varname = 'zthe'
417      zthe(:) = 0.0
418      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zthe,0.0,
419     .     jjm, rlonu, rlatv , interbar )
420      varname = 'zpic'
421      zpic(:) = 0.0
422      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zpic,0.0,
423     .     jjm, rlonu, rlatv , interbar )
424      varname = 'zval'
425      zval(:) = 0.0
426      CALL startget(varname,iip1,jjp1,rlonv,rlatu,klon,zval,0.0,
427     .     jjm, rlonu, rlatv , interbar )
428c
429      rugsrel(:) = 0.0
430      IF(ok_orodr)  THEN
431        DO i = 1, iip1* jjp1
432         rugsrel(i) = MAX( 1.e-05, zstd(i)* zsig(i) /2. )
433        ENDDO
434      ENDIF
435
436
437C
438C lecture du fichier glace de terre pour fixer la fraction de terre
439C et de glace de terre
440C
441      CALL flininfo("landiceref.nc", iml_lic, jml_lic,llm_tmp, ttm_tmp
442     $    , fid)
443      ALLOCATE(lat_lic(iml_lic, jml_lic), stat=iret)
444      ALLOCATE(lon_lic(iml_lic, jml_lic), stat=iret)
445      ALLOCATE(dlon_lic(iml_lic), stat=iret)
446      ALLOCATE(dlat_lic(jml_lic), stat=iret)
447      ALLOCATE(fraclic(iml_lic, jml_lic), stat=iret)
448      CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp
449     $    , lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid)
450      CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp
451     $    , 1, 1, fraclic)
452      CALL flinclo(fid)
453C
454C interpolation sur la grille T du modele
455C
456      WRITE(*,*) 'dimensions de landice iml_lic, jml_lic : ',
457     $    iml_lic, jml_lic
458c
459C sil les coordonnees sont en degres, on les transforme
460C
461      IF( MAXVAL( lon_lic(:,:) ) .GT. 2.0 * asin(1.0) )  THEN
462          lon_lic(:,:) = lon_lic(:,:) * 2.0* ASIN(1.0) / 180.
463      ENDIF
464      IF( maxval( lat_lic(:,:) ) .GT. 2.0 * asin(1.0)) THEN
465          lat_lic(:,:) = lat_lic(:,:) * 2.0 * asin(1.0) / 180.
466      ENDIF
467
468      dlon_lic(1 : iml_lic) = lon_lic(1 : iml_lic, 1)
469      dlat_lic(1 : jml_lic) = lat_lic(1 , 1 : jml_lic)
470C
471      CALL grille_m(iml_lic, jml_lic, dlon_lic, dlat_lic, fraclic
472     $    ,iim, jjp1,
473     $    rlonv, rlatu, flic_tmp(1 : iim, 1 : jjp1))
474cx$$$      flic_tmp(1 : iim, 1 : jjp1) = champint(1: iim, 1 : jjp1)
475      flic_tmp(iip1, 1 : jjp1) = flic_tmp(1 , 1 : jjp1)
476C
477C passage sur la grille physique
478C
479      CALL gr_dyn_fi(1, iip1, jjp1, klon, flic_tmp,
480     $    pctsrf(1:klon, is_lic))
481C adequation avec le maque terre/mer
482c      zmasq(157) = 0.
483      WHERE (pctsrf(1 : klon, is_lic) .LT. EPSFRA )
484          pctsrf(1 : klon, is_lic) = 0.
485      END WHERE
486      WHERE (zmasq( 1 : klon) .LT. EPSFRA)
487          pctsrf(1 : klon, is_lic) = 0.
488      END WHERE
489      pctsrf(1 : klon, is_ter) = zmasq(1 : klon)
490      DO ji = 1, klon
491        IF (zmasq(ji) .GT. EPSFRA) THEN
492            IF ( pctsrf(ji, is_lic) .GE. zmasq(ji)) THEN
493                pctsrf(ji, is_lic) = zmasq(ji)
494                pctsrf(ji, is_ter) = 0.
495            ELSE
496                pctsrf(ji,is_ter) = zmasq(ji) - pctsrf(ji, is_lic)
497                IF (pctsrf(ji,is_ter) .LT. EPSFRA) THEN
498                    pctsrf(ji,is_ter) = 0.
499                    pctsrf(ji, is_lic) = zmasq(ji)
500                ENDIF
501            ENDIF
502        ENDIF
503      END DO
504C
505C sous surface ocean et glace de mer (pour demarrer on met glace de mer a 0)
506C
507      pctsrf(1 : klon, is_oce) = (1. - zmasq(1 : klon))
508
509
510      WHERE (pctsrf(1 : klon, is_oce) .LT. EPSFRA)
511          pctsrf(1 : klon, is_oce) = 0.
512      END WHERE
513
514      if (couple) pctsrf(1 : klon, is_oce) = ocemask_fi(1 : klon)
515
516      isst = 0
517      where (pctsrf(2:klon-1,is_oce) >0.) isst = 1
518C
519C verif que somme des sous surface = 1
520C
521      ji=count( (abs( sum(pctsrf(1 : klon, 1 : nbsrf),dim=2))-1.0)
522     $    .GT. EPSFRA)
523      IF (ji .NE. 0) THEN
524          WRITE(*,*) 'pb repartition sous maille pour ',ji,' points'
525      ENDIF
526
527!      where (pctsrf(1:klon, is_ter) >= .5)
528!        pctsrf(1:klon, is_ter) = 1.
529!        pctsrf(1:klon, is_oce) = 0.
530!        pctsrf(1:klon, is_sic) = 0.
531!        pctsrf(1:klon, is_lic) = 0.
532!        zmasq = 1.
533!      endwhere
534!      where (pctsrf(1:klon, is_lic) >= .5)
535!        pctsrf(1:klon, is_ter) = 0.
536!        pctsrf(1:klon, is_oce) = 0.
537!        pctsrf(1:klon, is_sic) = 0.
538!        pctsrf(1:klon, is_lic) = 1.
539!        zmasq = 1.
540!      endwhere
541!      where (pctsrf(1:klon, is_oce) >= .5)
542!        pctsrf(1:klon, is_ter) = 0.
543!        pctsrf(1:klon, is_oce) = 1.
544!        pctsrf(1:klon, is_sic) = 0.
545!        pctsrf(1:klon, is_lic) = 0.
546!        zmasq = 0.
547!      endwhere
548!      where (pctsrf(1:klon, is_sic) >= .5)
549!        pctsrf(1:klon, is_ter) = 0.
550!        pctsrf(1:klon, is_oce) = 0.
551!        pctsrf(1:klon, is_sic) = 1.
552!        pctsrf(1:klon, is_lic) = 0.
553!        zmasq = 0.
554!      endwhere
555!      call gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque)
556C
557C verif que somme des sous surface = 1
558C
559!      ji=count( (abs( sum(pctsrf(1 : klon, 1 : nbsrf), dim = 2)) - 1.0 )
560!     $    .GT. EPSFRA)
561!      IF (ji .NE. 0) THEN
562!          WRITE(*,*) 'pb repartition sous maille pour ',ji,' points'
563!     ENDIF
564
565      CALL gr_fi_ecrit(1,klon,iim,jjp1,zmasq,zx_tmp_2d)
566      write(*,*)'zmasq = '
567      write(*,'(96i1)')nint(zx_tmp_2d)
568      call gr_fi_dyn(1, klon, iip1, jjp1, zmasq, masque)
569      WRITE(*,*) 'MASQUE construit : Masque'
570      WRITE(*,'(97I1)') nINT(masque(:,:))
571
572
573
574C Calcul intermediaire
575c
576      CALL massdair( p3d, masse  )
577c
578
579      print *,' ALPHAX ',alphax
580
581      DO  l = 1, llm
582        DO  i    = 1, iim
583          xppn(i) = aire( i, 1   ) * masse(  i     ,  1   , l )
584          xpps(i) = aire( i,jjp1 ) * masse(  i     , jjp1 , l )
585        ENDDO
586          xpn      = SUM(xppn)/apoln
587          xps      = SUM(xpps)/apols
588        DO i   = 1, iip1
589          masse(   i   ,   1     ,  l )   = xpn
590          masse(   i   ,   jjp1  ,  l )   = xps
591        ENDDO
592      ENDDO
593      q3d(iip1,:,:,:) = q3d(1,:,:,:)
594      phis(iip1,:) = phis(1,:)
595
596C init pour traceurs
597      call iniadvtrac(nq)
598C Ecriture
599      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
600     *                tetagdiv, tetagrot , tetatemp              )
601      print*,'sortie inidissip'
602      itau = 0
603      itau_dyn = 0
604      itau_phy = 0
605      iday = dayref +itau/day_step
606      time = FLOAT(itau-(iday-dayref)*day_step)/day_step
607c     
608      IF(time.GT.1)  THEN
609       time = time - 1
610       iday = iday + 1
611      ENDIF
612      day_ref = dayref
613      annee_ref = anneeref
614
615      CALL geopot  ( ip1jmp1, tpot  , pk , pks,  phis  , phi   )
616      print*,'sortie geopot'
617     
618      CALL caldyn0 ( itau,uvent,vvent,tpot,psol,masse,pk,phis ,
619     *                phi,w, pbaru,pbarv,time+iday-dayref   )
620       print*,'sortie caldyn0'     
621      CALL dynredem0("start.nc",dayref,phis,nqmx)
622      print*,'sortie dynredem0'
623      CALL dynredem1("start.nc",0.0,vvent,uvent,tpot,q3d,nqmx,masse ,
624     .                            psol)
625      print*,'sortie dynredem1'
626C
627C Ecriture etat initial physique
628C
629      write(*,*)'phystep ',dtvr,iphysiq,nbapp_rad
630      phystep   = dtvr * FLOAT(iphysiq)
631      radpas    = NINT (86400./phystep/ FLOAT(nbapp_rad) )
632      write(*,*)'phystep =', phystep, radpas
633cIM : lecture de co2_ppm & solaire ds physiq.def
634c     co2_ppm   = 348.0
635c     solaire   = 1365.0
636
637c
638c Initialisation
639c tsol, qsol, sn,albe, evap,tsoil,rain_fall, snow_fall,solsw, sollw,frugs
640c
641      tsolsrf(:,is_ter) = tsol
642      tsolsrf(:,is_lic) = tsol
643      tsolsrf(:,is_oce) = tsol
644      tsolsrf(:,is_sic) = tsol
645      snsrf(:,is_ter) = sn
646      snsrf(:,is_lic) = sn
647      snsrf(:,is_oce) = sn
648      snsrf(:,is_sic) = sn
649      albe(:,is_ter) = 0.08
650      albe(:,is_lic) = 0.6
651      albe(:,is_oce) = 0.5
652      albe(:,is_sic) = 0.6
653      alblw = albe
654      evap(:,:) = 0.
655      qsolsrf(:,is_ter) = 150
656      qsolsrf(:,is_lic) = 150
657      qsolsrf(:,is_oce) = 150.
658      qsolsrf(:,is_sic) = 150.
659      do i = 1, nbsrf
660        do j = 1, nsoilmx
661          tsoil(:,j,i) = tsol
662        enddo
663      enddo
664      rain_fall = 0.; snow_fall = 0.
665      solsw = 165.
666      sollw = -53.
667      t_ancien = 273.15
668      q_ancien = 0.
669      agesno = 0.
670cIM "slab" ocean
671      tslab(1:klon) = tsolsrf(1:klon,is_oce)
672      seaice = 0. !kg m-2
673c
674      frugs(1:klon,is_oce) = rugmer(1:klon)
675      frugs(1:klon,is_ter) = MAX(1.0e-05, zstd(1:klon)*zsig(1:klon)/2.0)
676      frugs(1:klon,is_lic) = MAX(1.0e-05, zstd(1:klon)*zsig(1:klon)/2.0)
677      frugs(1:klon,is_sic) = 0.001
678      fder = 0.0
679      clwcon = 0.0
680      rnebcon = 0.0
681      ratqs = 0.0
682      run_off_lic_0 = 0.0
683
684c
685c Avant l'appel a phyredem, on initialize les modules de surface
686c avec les valeurs qui vont etre ecrit dans startphy.nc
687c
688      dummy(:,:) = 1.0
689      ocean = "slab"
690      pbl_tke(:,:,:) = 1.e-8
691      zmax0(:) = 40.
692      f0(:) = 1.e-5
693      ema_work1(:,:) = 0.
694      ema_work2(:,:) = 0.
695
696      call fonte_neige_init(run_off_lic_0)
697      call pbl_surface_init(qsol, fder, snsrf, qsolsrf,
698     $     evap, frugs, agesno, tsoil)
699      call ocean_slab_init(dummy(1,1), tslab, seaice, dummy(:,:))
700
701      call phyredem("startphy.nc",phystep,radpas,ocean,
702     $    latfi, lonfi, pctsrf, tsolsrf,
703     $    albe, alblw, rain_fall, snow_fall, solsw, sollw,
704     $    radsol,
705     $    zmea, zstd, zsig, zgam, zthe, zpic, zval, rugsrel,
706     $    t_ancien, q_ancien, rnebcon, ratqs, clwcon, pbl_tke,
707     $    zmax0, f0, ema_work1, ema_work2 )
708
709
710
711C     Sortie Visu pour les champs dynamiques
712      if (1.eq.0 ) then
713      print*,'sortie visu'
714      time_step = 1.
715      t_ops = 2.
716      t_wrt = 2.
717      itau = 2.
718      visu_file='Etat0_visu.nc'
719      CALL initdynav(visu_file,dayref,anneeref,time_step,
720     .              t_ops, t_wrt, nqmx, visuid)
721      CALL writedynav(visuid, nqmx, itau,vvent ,
722     .                uvent,tpot,pk,phi,q3d,masse,psol,phis)
723      else
724         print*,'CCCCCCCCCCCCCCCCCC REACTIVER SORTIE VISU DANS ETAT0'
725      endif
726      print*,'entree histclo'
727      CALL histclo
728      RETURN
729      !
730      END SUBROUTINE etat0_netcdf
731
Note: See TracBrowser for help on using the repository browser.