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

Last change on this file since 932 was 927, checked in by lmdzadmin, 17 years ago

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